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A General Method for the Computation of Cartesian Coordinates 
And Partial Derivatives of the Two-Body Problem 


A general solution for cartesian coordinates and partial derivatives 
of the two-body problem has been programmed as a double-precision 
Fortran 4 program for the IBM 7094. The compact subroutine provides 
an accurate and efficient computation of coordinates and partial derivatives 
for all cases of two-body motion. A derivation of all equations used by 
the subroutine is given for those interested in the formulation. A de¬ 
scription of the subroutine is also given for those interested in the 
details of the computation. 
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INTRODUCTION 


The classical two-body problem may be reduced to the solution 
of the differential equations 

dr JW 


dt 


dr 

dt 


= “M 


where is a known constant and 

r_ yj-f • 

in the magnitude^ of r. A solution of the two-body problem gives the 
position vector r^ and velocity vector i^at any time t in terms of the 
position vector r^ and velocity vector r^ at a reference time t • Many 
practical applications also require the partial derivatives of the com¬ 
ponents x, y, z of r and x, 
x - Yn' z n of ^ an <* x , y Q j i Q of r Q 

z and partial i vauvco . • • • 

3 (x, y, z, x, y, z) 

8 | ‘o'' , o' , o’VV i o l 

have the disadvantage that several different formulations must be avail¬ 
able and one must be selected depending on the particular values of x 


X '0- z o of r o. ln ? 4 

nates x, y, z, x, y 


y, z of r with respect to the components 
Zq of r^. Classical formulas for the coordi- 
and partial derivatives 


'o’ o’ V V 


Zq and n 


0’ 


General formulas are available for computing the coordinates x, y,z. 


x, y,z and partial derivatives 


9 (x, y, z, x, y, z) 

8( W VVW 


as well as 


(x,y, z, x, y, z) , ,, , , , ... 

^ ^ - - for all possible values of x^y , z Q , x , y 0> z t , t and n . 

The formulas have been programmed as a double precision Fortran 4 sub¬ 
routine for the IBM 7094. The subroutine is advantageous in that it offers 
an efficient and compact program for the accurate computation of coordinates 
and partial derivatives for all possible cases of two-body motion. The 
purpose of this report is to give a derivation of all equations used by the 
subroutine and a detailed description of the computations. For those who 
are not interested in these details but only in using the subroutine, the 
inputs and outputs are described in Paragraph 4. 1 below and a Fortran 4 
listing is given in Paragraph 4. 8. 


The solution for the coordinates is derived in Section 1 by defining a 
new independent variable ip by the differential equation 


d ip 

dt 


r 
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which Sundman applied to a theoretical investigation of the three-body 
problem. The solution is a modification of that given originally by 
Stumpff and differs primarily in that it is valid for all values of the 
constant n . The relation of the general solution for coordinates to the 
classical formulas for elliptic two-body motion is also considered. It 
is shown that;/; is a generalization of the eccentric anomaly and that 
the relation between t and 4> is a generalization of Kepler's equation. 

The partial derivatives are obtained in Section Z by differentiating 
the solution for coordinates and manipulating the results algebraically. 
The approach is an extension of Sconzo's derivatives which are based 
on those of Kilhnert. Although the derivation is tedious and involved, 
the final derivatives are compact and general. Some formulas used in 
the derivation of the coordinates and partial derivatives are discussed in 
Section 3 with formulas that are used for computation. A detailed 
description of the double-precision Fortran 4 subroutine for the IBM 7094 
is given in Section 4. The description should be sufficient for those 
wishing to make modifications or additions to the subroutine. 


viii 



1* THE SOLUTION FOR THE COORDINATES 


1. 1 THE EQUATIONS OF MOTION OF THE TWO-BODY PROBLEM 

The equations of motion of the two-body problem are 

r — r 


m i r i 


- Gm i m 2 


1 


12 


rn r 
2 2 


- G m^m^ 


12 


where 

"12 ■ V <? > " ‘ ? 2> 

is the distance between the two bodies and G is the universal 
constant of gravitation. At the time t the body of mass m has the in¬ 
ertial cartesian position vector r^ and velocity vector ? , and 

the body of mass has the inertial cartesian position vector r 

- 2. 

and velocity vector 7 . Each body is attracted toward the other 

2 

with a force of magnitude G m^m^/r where the acceleration of 

m toward m is "r and the acceleration of m toward m is 7 . 

1 £ 1 2 12 

Initial conditions are specified by the position 7 and velocity 

* 1 ^ 

1,0 of m and the position 7 and velocity 7 of m at a 

2, U 2,02 

given reference time t . 

A closed -form solution of the differential equations above gives 

r , r and r , r at any time t in terms of r ,7 and 

1 1 . l c 1 , 01,0 

r 2 O' r 2 0 at t * me t Q as we ^ as the constants m^, m^ and G. Such 

a solution is extremely important, for many practical applications 

in astronomy and astronautics. However, it is convenient to trans- 

form the differential equations by defining the position R of the 

center of mass and the position 7 of m relative to m by 

2 1 


1 



m i r i + m z r z 

m l + m 2 


r = r — r 
2 1 


so that 


m i r i + m 2 r 2 
m i + m 2 


r 2 " r i 


give the velocity R of the center of mass and the velocity ? of 
relative to m^. Thus the values of the four quantities above 
at the reference time t^ are 

m, r* + rn rl ^ 

r = i jvQ 2 2,0 

0 m + 


V 


m, r, „ + m„ r_ 

1 1,0 _ 2 2 , 0 

m l + m 2 


r = r — r 

0 2,0 1,0 


r = r — r 

0 2,0 1 , 0 . 


Time differentiation of R and r above and substitution from the 
differential equations show that 


where 


"r - - G (m + m^) r/r~ 
r = V? • r 


An integration of these differential equations gives R, R and r, r 
at any time t. Simultaneous solution of the defining equations for 
R and r* shows that 


V R - 


_2_ 

m, + m. 


r 2 = R t 


m l +m 2 


2 



and time differentiation gives 


m„ 


R - 


m i 4 m 2 


m, 


r , = R + -——' 

£ m, + rn^ 


so that r* , r and r , r may always be determined from R, R 

11 L, Ld 

and r, r . 

Thus the integration of the differential equations for ? and r^ 
has been transformed to an integration of differential equations for 
R and Integration of the differential equation 

9l 

R = 0 

for R shows that 


and 


R = R„ + 


v- 


V 


so that the center of mass moves along a straight line at constant 

—» -V —► -V 

velocity. However, the inertial coordinates R^, R^ and thus R, R 
at any time t are usually unknown or the center of mass is chosen 
as the origin of inertial coordinates. In any case, the two-body 
problem reduces to the solution of the differential equation 

3 

r = - n r/r 

where 

tf=?—=*— 
r = v r • r 

is the magnitude of r and 
M = G (m^ + m^) 

is a postive constant. This differential equation is actually ob¬ 
tained for all inverse-square Newtonian forces between two bodies, 


3 



but the constant M may be different. For example, the same 
differential equation describes the electrostatic repulsion of 
two electrons but is then negative and is a more complicated 
function of the electronic mass and charge. 
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1.2 DERIVATIVES WITH RESPECT TO A NEW VARIABLE 


The differential equation 

- 3 

r = - M r / r 

where 

+ ^ 

r = r • r 

is conveniently integrated by defining a new independent variable 
4 by the differential equation of the Sundman transformation, 

$ = 1/ r 

where 4 J is zero when t is t^. Using a prime to denote differentiation 
with respect to 4> > 

t' = 1/4 = 1 / ( 1 / r) 


from the properties of derivatives. The chain rule for differentiation 

shows that 

/ -L / 
r = r t 

= r r. 

The method of solution described below requires determination 
of the successive derivatives of t and r with respect to 4 J . The 
second derivative of t is 


1 


1 


r' = — (r * r) 2 2r 


r r 


where the quantity 
a = "T • ■? 

is defined so that 

t" = r' = (J . 

The second derivative of r^is 

-V / -V / 

r = r r + r r 
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but 


r 

r 


= 't t' = (- nr/r 3 ) r 
2 

= - M r/r 


so that 


"r* = - /x r / r + o’ r. 
The third derivative of t is 


. 

t = r 


ff = r • r + r * r 


= rr • r + r • (- n r7 r ) 
-V ^ 

= r (r • r) - \i . 


However, the quantity 
# » 

O' = r * r - 2 fx / r 
is defined to obtain 

t"' S r " - a= r(r -7) -2 n + fi = r(r • "r - 2 /x / r) + n 

= a r + n . 

The third derivative of r is then 

-^/ -► / . 2 /*^ , / 
r =-/xr/r+/irr/r + (Jr + ar 

-V -► 2 -U. -► 2 

= - fi r + /i r a/r + (ar+^)r-a,^r/r 

= o r r. 

The quanity a has the important property that 

9 ^ 2 

o; / = 2r • r' +2 (ir'/r 

A- -► 2 2 

= 2r -(-/ir/r )+ 2 |i a/r 

2 2 
= -2 \i a/r + 2 a/r 

= 0 


so that cl is constant for every value of ip as well as t. This 
property gives all the derivatives of t and r as simple functions 
of the first two derivatives t 1 , t" and r 1 , r n . From the results 
above, 

t m = a t* + m 
r ,tt = a r* 1 


6 



so that 


t"" = at" 

J*iui _ Q ~£<i 

t’"" = a t'" = a ( a t' + n ) 

2 , 

= a t' + a n 

r ,im = a ~r'" = a (a f 1 ) 

2 

= a r 1 

t" H " = a 2 t" 

—► 2 -*■ 

jm HUM _ Q| r 1l 

t"""' = a 2 t'" = a Z ( a t' + n ) 

3 , 2 

= a t' + or n 

J.IIIIMI - Qr ^111 - Of ( Of 7') 

3 - 

= a r ’ 

etc. 

The final result of the above derivation is that 

t' = r 

t" = o 

t m = a r + n 

t ,,n = a (J 

2 

t ,,m = a r + a n 

t"""= a Z a 

t ,, ""i = a Z r + a fx 

t lilt'll! = (% ^ a 

etc. 


7 



and that 


r' 

= 

♦ 

r ~r 

7«, 

= 

- n r/r + CT r - 

—► 

j. Ml 

= 

a ( r r) 



a (- fz "r / r + a "r) 

-jflltll 

= 

2 . ^ 

Q? (r r) 

r umi 

= 

(_ /i 7/r + o’~r) 

r mini 

= 

3 -i. 

a (r r) 

■^11111111 

= 

3 _ 

a (-jur/r + ar) 


etc. 


where 


a = r. 


-2 n /: 


since a is zero and a is invariant for all values of # as well as t. 
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1.3 SERIES SOLUTION IN THE NEW VARIABLE 


The Taylor's series for t in terms of ip is 

t = l o + V * + t o" ^ 2/2: + V" ^ 3/3! + V" ^ 4/4! _ 

where t Q , t Q ', tg", t 0 "', t "". are the values of t, t', t", 

t" 1 , t"", ... when ip is zero, i. e. when t is t . But 


t o" 


= a 


V" 


■ “ r 0 + <* 


t lilt 
0 


t mn 
0 


= Oi (T 


0 

2 

a r Q + q t ix 


t Mmi = a a 
0 0 


etc. 


from the results above where 

“ = W 2 " /r o. 

Thus substitution for t^', t^", t^'”, t^ ,m , . 
series gives 


in the Taylor’s 


* - *o + r o * 


+ «■„ * /2! 


3 4 

+ (o'r 0 +fi)i^/3!+o; tp / 4! 

2 5 2 6 

+ ( a r^-tan) tp / 5! + Of cr^ ip / 6! 
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or, collecting coefficients of r^, o^, and pi , 

t = t 


3 2 5 3 7 

+ r^(^+a^/3!+o! ip / 5 ! + a ip I 7 ! + ....) 


+ <r^ ( ^ 2 /2! + a ip 4 ! + O' 2 tf^/6 ! + a? 3 0 8 / 8! + 


+ pi ( ^ ^/3 ! + Of i/ ^/5 ! + if 2 


7 3 9 

i/> / 7 ! + a ^ /9! + 


gives the value of t corresponding to any particular value of 
Similarly, the Taylor's series for r in terms of ip is 


where 


' = ? o + ? o* + "o' * 2/2! + " 0 ” 


* /3! + r Q "" 


ip /4 ! + . . . 


V 


= r r 
0 0 


- “ V r o * °o *0 


= ct r r 
0 0 


a (- n r 0 /r Q + r ft ) 


0 O' 


_2 - 
a r r 
0 0 


f 0. * “ <- 'V'o* "oV 

etc. 

from the results above. Thus substitution for “r^ 1 , r^", r^ M \ 
in the Taylor's series gives 


? - ? o + r o r o * 


+ <* V r o + ’o'o 1 * /2! 


+ a r Q r/3! + « (- M^/r^ ^ T) * /4! 

+ a 2 r 0 Tq /5 ! + « 2 (- M r Q + °q r Q^ 


or, collecting coefficients of and r^, 
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r =[1 - 


( 4) Z l2l + a ip*/ 4 ! + ot ip u /6\ + a ip"/ 8 ! + 


)lr, 


0 


r ( ip + a i/i 3 /3 ! + a ip ^ / 5} + a'* ^ 7 / 7 ! + ...) 

0 

+ cr 0 ( 0 2 /2! + a ip 4 / 4! + a 2 ip b / 6! + a 3 ^ 8 / 8! + ...) 

which gives the value of "r corresponding to any value of # . 

The results above are conveniently summarized by defining the 
transcendental functions 


S 2 = 


S 3 = 


1 + 

a 


/ 2 ! 

+ 

2 

a 


4 / 4! 

+ 

3 

a 

ip 

6 / 6! 

+ . ., 

Ip + 

a 

>P 3 

/ 3! 

+ 

2 

a 

>p 

5 / 5! 

+ 

a 3 

ip 

7 / 7 ! 


Ip 

Z /2! 

! + 

a 

* 

4 / 4! 

+ 

2 

a 

>P 

6 /6! 

+ 

3 

a 

ip 8 / 8! 

Ip 

3 /3 

! + 

a 

ip 

5 /5! 

+ 

a 2 


? /7! 

+ 

3 

a 

* 9 /9! 


which have the derivatives 


V = as i 


2 1 

V = S 2 

with respect to ip . In this notation 

* - *0 + r 0 S l + a o’z* “ S 3 

gives t" for any value of ip . Also 

r = f = r Q s j 1 + a Q s 2 ‘ + \i s^ 

= rs+ O’ s, + H s 
0 0 0 1 2 


and 


r ' - r 0 V + % V + " S 2 ,= r 0 “ S l + "o*0 + " S 1 


= Vo + 1 “ r o + M)S 1 
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gives r and cr corresponding to any value of if). 
It is also convenient to define the functions 


f = 1 


M S 2 /r 0 


* ' r 0 =1 + °0 S 2 
which have the derivatives 

f = - « s 2 V r 0 
= f V r o 

8' = r 0 S l' + Vi' 


= r s„ + cr s , 

0 C 0 1 


However, the equation for t as a function of ip shows that 


1 - 4 0 - M S 3 


is an alternate expression for g. Differentiation of this equation 


shows that 


g 1 = t 1 - ix 
= r - V S o 


is an alternate expression for g 1 . 

In terms of f and g, 

r = f ~r + g r 
0*0 

gives 7 corresponding to any value of if). Differentiation of this 

equation with respect to ip gives 

r = f r + v r 
0 * 0 

and differentiation of "r with respect to t gives 

7 = f r + g r* 

0*0 


12 



where 


and 


f = f <P = (- M Sj/Tq) (l/r) 

= -C W> 


g = g' $ = (r - s ) (1/r) 

= 1 - /i s 2 /r 


or 


g = g’ tp = (r s + cr s) (1/r) 
s 5 r x 0 0 0 1 

= (r s)/r, 

o o or 

Another form of the equation relating t and ip above is obtained 
by defining 

y = r (r * r) - fi 
which satisfies the identity 

y = r (? • 7 - 2^i / r) + n 


Then 


= a r + n . 


1 ■ V r 0 S 1 + % S 2 + “ S 3 


so that 


■ V r 0 < * * “ s 3 )+ Vz + M8 3 

' V r 0 * + % s 2 + < “ r 0 + " ) s 3 


* * Vo * + < ’0 S 2 + T 0 S 3- 


Also, differentiation of this equation gives alternate expressions 
for r and o . 


/ 

t = r 


r + cr s, + y 
0 0 1 r 0 2 


% s 0+ y o S] . 
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In addition 


-/ -V 


u/ 


V' = (r . r)' = r . r + f -r = r r . r + 7-(- „ r/r“) 

= y 


so that 


t in = r 1 ’ - a ' ~ y = y s + cr or s 

' r 0 0 0 1 


= r(t 
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1.4 GENERAL SOLUTION FOR COORDINATES 

The formulas above may be used to calculate r, r from given 
values of "r , r^ and \i , t Q> t. First 


VT • ?„ 
0 0 


r = r • r 
0 0 0 

“ = V'o - 2 " /r 0 


are determined. Then the parameter ip and its transcendental functions 


s o 

1 + a 

2 2 4 3 

ip /2! + a ip /4! + a 

ip^ / b\ + .... 

S 1 = 

S 2 = 

ip + a 

3 2 5 3 

ip /3!+ a ip / 5! + a 

j/> 7 / 7 ! + . . . . 

P Z !Z\ 

+ a ip 4 / 4! + a 2 ^/6! + 

a 3 i/) 8 /8 ! + . . . 

s * = 

^ 3 / 3 ! 

+ a ^^5! + a 2 ip7/l\ + 

3 9 

a $ / 9 ! + .. .. 


are obtained by solving the equation 


for tp 


and 


0 

Then 


t + r s + <r s 0 + ns 


0 1 


0 2 


r = Vo + 5 1 4 “ S 2 


1 - “ s z /r o 


g 


(t - t 0 ) - X s 3 


f = - n s l /(r r Q ) g = 1 - M s 2 / r 

give the final solution for the coordinates. 

' = ‘ T o 4 8 7 0 

? = i ?o + g 

These are the equations for coordinates given previously by the 
author (1965). 
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The great advantage of these formulas is their complete 
generality. They are a modification of equations which were 
derived by Karl Stumpff (1947). Stumpff's derivation differs 
somewhat from that above and he normalizes his equations for 
computational convenience and simplicity. Similar formulas 
have been obtained independently by Herrick (I960) and Sperling 
(1961), but Stumpff was apparently the first to recognize that a 
generalized solution for coordinates of the two-body problem 
could be obtained by utilizing new transcendental functions sim¬ 
ilar to Sq, s i , s^, s^ above. Stumpff 1 s work (1947) is given in 
his textbook (1959) and is also available in English (1962). The 
formulas above differ primarily in that the equations are not 
normalized but are left in a form that is continuous through 
the trival cases where n or (t - t ) is zero. The formulas are 
thus equally valid for negative values of M or (t - t ). 

The transformation from the variable t to the variable ip 
is often treated as a regularizing transformation since the sol¬ 
ution for 7 above is continuous through a collision which can 

occur if n > 0 and r is parallel to . For such a collision 

U 0 

r is indeterminate at the origin but ?* is zero at the origin and 
continuous through the collision. Also, r and a equal zero at 
the origin and are continuous through the collision. The con¬ 
tinuation of the solution through a collision is of course of 
little physical importance. However, it is practically advantageous 
because numerical problems are eliminated in the computation of 
position coordinates for near-collisions. 

The regularizing transformation from t to i p was used by 
Sundman (1912) in an investigation of the three-body problem. 
Stumpff realized that the transformation was of computational 
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value in solving the two-body problem, and that the resulting 
formulas could be placed in a form that was equally valid for 
elliptic, parabolic, and hyperbolic orbits including the cases 
of circular and rectilinear motion. He also realized the re¬ 
lation of the general formulas to classic formulas for elliptic, 
parabolic, and hyperbolic motion. This provides an alternate 
way of deriving the general solution, i.e. it may be obtained 
by manipulating classical formulas. Herrick (I960) used this 
approach and his formulas were modified to obtain the general 
solution above before Sconzo informed the author of Stumpff's 
prior work. 
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1.5 CORRESPONDENCE TO CLASSICAL FORMULAS 
In the case where a is negative, let 

-JL- = _ 

a 

define the length a and let 

AE = M /a i p 

define the angle AE. Then the functions s^, s^, s^, may be 
expressed in terms of y/Wa and trigonometric functions of AE. 

For example, 

s 3 = ip 3 / 3! + a 4> S /5\ + a z ^ 7 /7! + a 3 ^ / 9!+.... 

= ip 3 / 3! - ( n /a) ^ 5 / 5! + ( H /a) 2 ^ ? /7! - ( M /a) 3 ^ 9 /9! + 
=y /^/a 3 i/) 3 /3!-V f q/a 3 $ /51+y/q/a /7! .... 

\fpja ~~ 3 

= (\ryjTjP ) 3 / 3 ! -(/77M0 5 /5! -1 ( \/ n/a >) 7 /7! - .... 

3 

y/ M /a ' 

= a 3 E/3! - A 5 E/5! + A 7 E/7 ! - .... _ 

3 

= AE - ( AE - A 3 E/3! + A 5 E/5! - A^E/7! - -) 

/ M/a 3 

= A E - sin A E 
v^TT- 3 

with similar derivations for s^, s^, s^. 
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The result is that 


s Q = cos AE 

sin AE 


1 Jjr> 


1 — COS AE 

( V /a) 


A E - sin AE 


S 3 = 


3 

These formulas give the following alternate solution when a 
is negative (which can occur only when p is positive). First 

r o 1 y ? o • r o 


a = 
0 


r * r 
0 0 


( p/a) = 2 n lr Q - r 0 * r Q 
are determined. Then solution of the equation 


t = t rt + r 


sin A E 


o o / 7 

\f M/< 

gives A E. Then 


r = r^ cos A E -I* 


+ a 0 


1 - cos AE 

( /*/a) 


+ M 


AE — sin A E 


sin A E 


° /i 


and 


f = 1 - 


f = - — 

r r 


1 - cos AE 
/ a 

sin A E 


give 


o v/m /a 
r = f "? A + g r 




+1* 


8 ' ‘ " ‘o - " 


1 — cos AE 

( M / a ) 

AE - sin AE 


= 1 - 


fx 1 — cos A E 


( n /a) 


r = f r 0 + » r 0 
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This solution is also valid when a is positive although AE and 
y^/ju /a are then imaginary quantities. However, the solution is 
not valid when a is zero and its numerical accuracy decreases 
as Oi approaches zero. 

The formulas for t, r, f, g, t, g above all reduce to classical 
forms such as those given by Herget (1948) if 

ae = e-e 0 

where E and E^ are, respectively, the eccentric anomalies at 
time t and at time t^. The well-known relations 

r Q = a < 1 - ecos E ) 


% 


a e sin E„ 


0 


are used to obtain the results.For example, 

sin AE 1 - cos AE 


t = t + r 
0 0 


4 o 


VV/i 


0 


M /a 




AE - sin AE 




t^ + a (1 - e cos E^) sin AE / t/ pi /a 


,j^a e sin (1 - cos AE) / ( p /a) 


+ n (AE - sin A E) / \J /a 
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so that 



(t 


t Q ) - (1 - e cos Eg) sin AE 

+ e sin Eq (1 — cos AE) 

+ AE - sin AE 

= sin AE - e cos E^ sin A E + e sin E 

0 0 

- e sin Eg cos Ae 


+ AE - sin AE 


= AE + e sin Eg — 


= AE + e sin E 

0 


e (sin Eg cosAE + cos E 
e sin (E Q + AE) 


= (E - E 0 ) + e sin E Q - e sin [ E Q + (E - E^) 

reduces to the classic form of Kepler's equation 

E - e sin E = E Q - e sin E Q + /„/a 3 (t - t Q ). 

Since all the steps above are reversible, an alternate derivation is 
to define 


, E - E 
0 = 0 


sjrr 

and express Kepler's equation and the classic relations 
t = a (1 - e cos E) 


g = t - t 0 - 


1 - - [ 1 - cos (E - EJ | 

r o ® 

(E - E o ) - sin (E E q ) 




f = - 




r o 


g = 1 - 


a 

r 


sin (E - E ) 

0 


[ 1 - cos (E - E q )] 


sin AE) 


1 
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in terms of ip and a . 

In the case where a is positive, let 

AF = \Z~a~ 4> 

define the angle A F. An equivalent definition 
AF = i AE 
where i = \Ti and 

a = - M / a * 

Either definition gives 

s. 


0 


cosh A E 
sinh F 


v a 

cosh AF — 1 

a 

sinh AF - AF 


s 3 ' 


\/ ~a 


T 


which may be reduced as above to obtain classic formulas for 
hyperbolic orbits. The solutions in terms of A E and AF are 
both equally valid for all cases where ot is not zero, but im¬ 
aginary quantities are eliminated by utilizing AE for negative 
a and AF for positive a . However, both solutions decrease 
in accuracy as Ot approaches zero. 

In the case where a is zero, the functions s^, s^, s^, s^ 
reduce to 

s o 1 1 

s l = $ 
s 2 = ip 2 /2 
s 3 = i/) 3 /3 ! 
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since only the leading term in these series is then non-zero. The 
resulting formulas may also be reduced to classical formulas for 
parabolic orbits. This suggests the possibility of using AE for 
negative a , ip for zero a , and A F for positive a , which 
corresponds to the classical approach. However, use of ip in the 
general solution for all values of cl is advantageous for several 
reasons. Computer storage requirements are reduced since only 
one formulation is required. The general soltuion is also contin¬ 
uous through zero ^ so that small values of a do no require a 
separate formulation for nearly-parabolic orbits. These advantages 
can become very important in practical applications where diff¬ 
erent types of orbits may be encountered. 
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2. THE PARTIAL DERIVATIVES OF THE SOLUTION 


2. 1 DIFFERENTIATION OF THE SOLUTION FOR COORDINATES 

The solutions r and ~r are functions of the initial conditions "r 

0 

and as well as p and the times t and t^. The differential re¬ 
lationships giving dr" and dr in terms of dr^ and dr^ as well as dpi 
and the time differentials dt and dt^ are very important. The 
differential relation between dF and dr^, d?" is given by Herget in 
a form due to Bower (1932) which is similar to the earlier work 
of Kilhnert (1879). Sconzo (1963) has extended this approach to give 
the differential relation between dr and dr^, dr^ as well. The 
approach given here is similar but concise expressions are ob¬ 
tained which are valid for all cases of two-body motion. The 
differential relationships are obtained by first differentiating the 
equations of the solution for coordinates and then combining the re¬ 
sults to eliminate all differentials other than dr", dr*" and dT^, dr" , 
dfi, dt, dtg. 

The differentiations are as follows. 

Differentiation of 


gives 


-*• r -* . -*■ 

r = f r 4 g r 

0 5 0 

-!► • - .4 

r = f r + g r 

0 5 0 

d? = f dr 0 + g d"? Q 

+ 7 o d ‘ * T o d 8 

= * d7 0 + 8 d7 o 

+ ? 0 d ' + T 0 d 8 
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Then 


gives 


and 


gives 


f = 1 - ^s i /r Q 


<*£=-'„ s 2 /r 0 - fi d. 2 /r 0+ „s 2 d r 0 /r 0 ‘ 

= [-(£ - 1) d r Q - „ d s 2 - s 2 d„|/ r 0 


g - t - t 0 - d = 3 


dg = dt - d t Q — d n s 3 - |i d s 3 


But 

= -/ids 3 -s 3 d M 
the alternate expression 

+ dt - dt^. 



g = r s + 

6 0 1 

cr s 

0 2 



for 

g gives 





dg = r d s 
b 0 1 

+ s i d r o + 

d 0 d s 2 + s 2 

da o 


= s i d r o 

+ r o ds i + 

"0 d 3 2 + =2 

d a o 


Similarly, 


gives 


‘ - - “ a A r r o) 

df = - d„ .,/frrJ- „ dSj/^r^t „ 8j dr 0 {r r^+KS, d 

= - ‘ dr 0 /r 0 “ ‘ dr/r d S l^ r 0 )- d “ V( r 'o') 


and 


g = 1 -fi s 2 ^ T 

gives 

i 2 

dg = - d n s /r - n d s /r + ^s dr/r 

^ L, L, 

- [-(id - (g- 1) dr - S 2 djU ^ ^ r - 
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But the alternate expression 


* = (r o s o + °0 5 i )/r 


gives 


d£ - (r 0 ds Q + s Q dr 0 + dSj + Sj d 


- ,r o s o* < ’o s i )dr/r 

■ l s 0 dr o - * dr * r o ds o * **o d9 1 * s 1 dd o l/r - 

The differentiation of the equation 
* * *0* *0 S 1 * % S 2* d S 3 

is equivalent to equating the two different expressions above for dg. 

This gives 

- „d S3 - s 3 d* + dt - dt p = .J dr 0+ r 0 d.j + o ds 2 * 3 2 d< 0 

r 0 d9 l * °0 ds 2 * <* ds 3 = - S 1 dr 0 - s 2 d % - S 3 d “ 

* dt - dt 0 

Similarly, differentiation of 

r = Vo* Vl* “ S 2 

is equivalent to equating the two expressions above for dg* This gives 
[- n ds 2 - (g-1) dr - s^d p] /r 

= [s o dr o - g dr + r 0 ds 0 + CT 0 ds i + s i d a o ] /x 

dr=5 0 dr 0* S l d< 0 + s 2 d 


* r 0 ds 0 * C 0 d9 l * “ dS 2. 
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The differentiation of 

s 0 = * + « # /2! + a ip 4 /4! + a 3 ip^/6! + .... 

S 1 = ^ + a !/i 3 /3!+ a 2 " ip $/ 5! + a 3 j^ 7 /7 ! + . . . . 

s 2 = # /2! + a i/) 4 /4!+Q! 2 ^ 6 /6!+ a 3 ^ 8 /8!+... 

s 3 = ^ 3 /3 ! + a i/' 5 /5!+Q! 2 i/) 7 /7!+a 3 ^ 9 /9! + .... 

is conveniently expressed by utilizing 

2 o / 

^/2!+ 2 Q/ ip / 4 ! + 3 Qf ip / 6! + . . . . 

ip / 3! + 2 a ip / 5 ! + 3 a 2 ^ 7 / 7 ! + .. .. 

^ 4 /4! + 2a ip ^ / 6\ + 3a 2 ^ 8 /8 ! + .... 

ip^ / 5\ + 2 a 7 / 7 ! + 3 a 2 ^ 9 /9!+.... 


3 s. 


9a 


9s 


1 


9a 
8s „ 


8a 


8s, 


8a 


Term by term differentiation of the series and combination of the 
results gives 


ds = 

0 

a s d f + 

0 

8a 

d a 

ds i = 

s 0 <**♦ • 

8s i 

8a 

d a? 

d ?2 = 

d * + - 

® S 2 

8a 

"da 

dS 3 = 

d ^ +- 

, dS 3 

8a 

-da 
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Finally, differentiation of 


r = (r • r ) 

0 ' 0 o' 


1/2 


CT 0 r Q ’ r 0 


“= V r o " 2 " /r o 


gives 


dr o = 


2 (? 0 • V 2 2? 0 d? 0 


= d ~f • / r 

0 0 0 


d O’ = r • dr + r - dr 
0 0 0 0 0 


da = 2 r Q • dr Q - 2d n /r Q + 2 n dr Q /r Q 

= 2 < Mr o /r o 2 * V d V d " /r o ) - 
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2.2 THE COMBINATION OF THE DIFFERENTIALS 

The combination of the differentials above to give dr and dr 
as a function of dr^, drd fx , dt, and dt^ is simplified by using 
the relations. 

? 0 = 8 7 - g * 

I f « 

? 0 = -f r I f r 

and 

fg - f g = i 

which are derived in Section 2.4 below. 

Substituting the two expressions for and r^ into the formulas 
above for d”? and dr* gives 

dr" = f Jr 4 g dr" 4 r" d f 4 r" dg 

= f dr^ + g + (g t - g r) df 

4 (-f "r 4 f *r) dg 

so that 

d? = f d? 0 + g dT 0 + r (g d f - 1 d g) 

+ r* (- gdf+fdg) 

and, similarly, 

d? = f d(? 4 g <fi* 4 r" d f 4 *?" d g 

0 & 0 0 0 s 

= f d? 0 + g d? 0 + (g r*- g r) d f 
+ (-f ? + f r) d g 
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so that 


dr = f dr" 0 + g dr^ + ~r (g d £ - fdg) 

+ r”(-g d f + fdg). 


The expressions (g d f - f d g), (-g d f 4 f d g) and 
(g d f - f d g), (-g d f 4 f d g) are evaluated in order as follows. 


Firstly, 

(gdf- fdg) = ?H f -1) dr Q - ^ ds 2 ' s 2 d M 1 / r 0 

- ila l dr 0 + r 0 dS l + % dS 2 + S 2 d< o' 


= [-f Sl -(f-1) g/r Q ] dr Q - f r Q d Sj -(g M /r Q + f o q ) 


ds. 


- fS 2 da 0 " ( * S 2 /r 0 ) ^ 


9 s 


1 


= l- f s i - g^o 1 dr 0 _fr 0 (s 0 d * +_ 53~ dQ; > 

9S 2 

-(g^ /r Q + f ff ) (s di^ + — da) 


" f s 2 d<T 0 “ S 2 /r 0 } d,i 


= [-fSj - (f- l)g/r 0 ] dr Q - [fr o s o + (g(i/y i ^JSjl 


dip 


9s. 


9 s. 


- [f r 


0 9o! 


+ (g M ! X n + f O' J 


9a 


da 


f s 2 d a 0 - (g s 2 /r 0 ) d M 


= M Sj - (f - 1) g / r Q ] dr Q - [f 


r s + a s 
0 0 0 1 


r r. 


g] r d ip 


Ms, 


- [- 


r 0 9a 


+ ( 


r s + cr s 
0 0 0 1 M 


1 * ”2 , 

a J —5^—] d a 


r i*o 0 9a 


- f s 2 d a Q - ( g s 2 /r Q ) 
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[ ~ f S 1 " (f - 1} g^o 1 dr 0 -(fg-fg)rd* 

8 s_ 


+ "r < s i 


9 S 

1 - s 

~3a~ 0 


~5cT > d « -f s 2 d a Q 


( g s 2 /r Q ) 


= r-f Sj -(f-1) g/r Q ] dr Q 


2(s 


_!!i 

1 da 


ds 


- s 


2_v do; 

c5a ' ? 


~ f s 2 da 0 - { § s 2 /r 0 ) 


But the identity 
2 

s 2 = 2 (s 


8s 


1 


3s. 


1 da 


- s 


0 Ba 


is treated in Section 3. 1. 

This gives 

<gdf - fdg) = [-f.j .{M| i/ Io , dr 


-Ms. 


2 ~ 


r S 2 2(Mdr o /r 0 + r 0 ' d ^0 " dM/r 0 )/2 
- f s 2 d ^ _(g s 2 /r o ) 

t-fs, -{£-l)g/r 0 Idr 0 - <g-l) * 2 («dr 0 /r 0 2 t T 0 .d? 0 

f S 2 d a 0 ~ ( g S 2 /r 0 ) ^ 

[-f Sj -(f-l)g/r Q + [{f-l) (g- 1)/ ro J dro 
-(g-1) s 2 • dr Q - f s 2 d <r Q - (s^i^) d/i 

[-£ D/r 0 ] d? 0 . r Q /r 0 - (g-1) ^ 7 . d ? o 

- £ s 2 (?-d? o + T g . d?) - (s 2 /r 0 ) dd 


d M / r Q ) 
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or, finally, 


(g d f - f < 

Secondly, 

(-g d f + f 


f s. + (f— 1)/ r 

i g) = - --- • dr"„ - f s, r*„ . dr - . 


0 


0 0 


2 0 0 


- d S 2 ? 0’ dr *0 - ( ®- U S 2 V d '0 " (S 2 /r 0> d " 


d g)= -g [-(f-l) dr Q - |jds 2 - s^d/i ] /r Q 


+f fs l dr 0 + r 0 dS l + 7 0 dS 2 +S 2 d V 


= If Sl + (f-1) g/r Q ]dr Q + f r Q ds } + (f O q + g Hfr Q ) ds 2 


+ f s 2 d a Q + (g s 2 /r 0 ) <** 


9 s, 


[f s i+ (f-1) g / r 0 ] dr Q + fr 0 (s 0 df + 

a 


da ) 


+ (f CT o + g ** / r Q ) dl/j + 


— da) + f s d°. 

da 2 0 


+(g s 2 /r Q ) 


t f S 1 + (f ^ g//r 0^ dr 0 + ^ r 0 S 0 + a Q S l) + gM s 1 / r Q l d ^ 


8s Bs 2 

+ tf qtt + ( f CT 0 + gfi /r n )_ ^r ] dQ; 


0 da 


+ f S dCT Q + (g S 2 /r Q ) d ^ 


[f Sj + (f-1) g/ r o 1 dr Q + (f 


r o s o + Vi 


r r 


1 g) r dip 


9 s i 


8 s - 


+ t* r n -^ +(f^ + gM/r 0 )-^l d« 


0 8a 


+ f s 2 d a Q + (g s 2 /r Q ) d M 
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= [f s j + (f — 1) g/r Q ]dr Q + (f g - f g) t dip 


9S 1 . 9s 2 

+ (f + (f CT 0 + g V-! r n ) — J d a + f s o d a n + (g s o /r n )dfi 


0 


da 


0 8a 


2 0 e 2 0' 


r 9s i 

= t fs, + a-1) g/r 0 l dr 0 + rd* + ^fr 0 — 

1 


+ (f a Q + gM /r o ) da H^d a Q + (g d M 


But 


r n ds , + dS o + M ds - 


0 ““1 • 0 '*“2 ■ ™3 " " s i dr o ~ s 2 d<7 0 ■ s 3 dM + dt ■ dt 0 

from above so that 

9S 1 9s 2 9s 3 

r 0 (s 0 d 4>+ ~W da > + % ( s i d «) +M < s 2 d ^ + da ) 


= - s l dr 0 - = 2 % " S 3 d " +dt - dt 0 

9s i 9s 2 9s 3 

s - (*« - + CT ~ - + **-— ) da 


0 0 0 1 2 


0 


da 


0 


da 


da 


- s. dr. - s d a- s d jx + dt — dt 


1 0 2 0 
9 s 


0 


r d ip = - (r 


0 9a 


1 8s 2 8s 3 

+ d o-sr +,, -^ ) d “ 


s l dr 0 ■ S 2 d °0 - S 3 d ‘‘ +d * - dt 0 


This gives 

(-g d f + f d g) 


9 s, 


= [ f Sj + (f- 1) g/r Q ] dr Q - (r Q 


9a 


+ % 


ds 9s 

— +M—) d a 

9a 9a 


- Sj dr Q - s 2 d a Q - s 3 dM+ dt - dt Q 


9s, 


+ t £ r o — + < f % * / r 0 > 

da 


9s, 


9a 


-] da 


+ f s 2 d % + ( g s 2 /r Q ) dM 
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[ (f-1) Sj + (f-1) g/r Q ] dr Q 


8 s 


- r. 


8a 


_ <r 


0 8a 
8 s 


8s 2 S , , 

- - n - + f r 


9a 


0 9a 


+ (£ V 8 " /r 0> 17 J 


dad (f-1) s 2 d <j Q + (g s 2 /r Q - s 3 ) djx 


+ dt - dt^ 

But the coefficient of da is 
8 s 


(f-1) r Q T-d- + [(f-D % + g^ 1 / r 0 ] 


d s. 


- M 


8a 


9 s, 


8a 


-M s. 


1 


v C. 


0 9a 
8 s 


+ [ 


-fJL S. 


0 


a + ( r s + a s ) M / r J 
0 0 1 0 2 0 Ba 


8 s 9 s, 

2 3 


1 


8a 


8S 2, . 8S 3 

- s -) + 

1 9a 9a 


IV] 

da J 


and the identity 


S 2 S 3 = 2 (s. 


d s 


da 


- 1 


9 s . 
s. 2 ) 


9a 


is treated in Section 3.1. This gives the following coefficient for da 

9 s 9 s 3 

-M < s 2 S 3 /2+ TT’ 0 (g S 2 ' 8 = 2 - “ S 2 S 3 " “ 2 H )/2 


8 s. 


= 16*2 - VoV 2 

8 s. 


8a 


-3 2 


s 2 - 2 9a 


9a y 


Thus 


s 

(-gdf+fdg) = [ (f-1) Sj +(f — 1) g / r 0 1 dr Q + (g 

r\ 

!) 


- 2 ir 


da / 2 + (f- 1) s 2 da Q 


da 


+ (g s 2 /r Q - s^) c%z + dt - dt Q 
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Setting 


8 


U = s 2 (t-t 0 )+ „.2 


da 


gives 


(-g d f + f d g) = [ (f- 1) Sj + (f- 1) g/r Q ] dr Q 

+ (g s 2 - U) ( M dr 0 /r Q 2 + * 0 . dr Q - d M/r Q ) 

+ U-l) s^d ct q + (g s^/i^ - s 3 ) dA* + dt - dt Q 
= [ (f“l) Sj + (f-l)g/r Q - (- j« s 2 /r 0 )g/r Q ] dr Q + U (-g/r Q 2 )dr o 


+ (g s - U ) r • dr + (f- 1) s d or 

B 2 0 0 2 0 


+ (-g s /r Q + u /r 0 + g s 2 ' r 0~ S 3* + at ~ dt 0 

= [ (f-l)s 1 + (f-l)g/r Q - (f-l)g/r Q l drj ) -‘r () /r 0 


+ U (- g /r 0 ) dr o T o /r o + (g s., -UJr^d^ 


+ (f-1) s 2 d cr Q + ( U /r Q - s 3 )d fi + dt-dt Q 


■ W V d? o + u ( "“ 7 o /r o >' d? o ~ u V d? o 


+ g s„ ~r • d7 + (f- 1) s fr • dr - + 7 • dr\) 
6 2 0 0 2 0 0 0 0 


+ ( U / r - s ) dg + dt - dt 
V ' 0 3 0 


and, finally 
(-g d f + f d g) = 


(f-1) s 


1 


r V dr 0 + (f ‘ 1) 3 2 r 0' dr 0 +S " 1|S 2 r 0' dr 0 


• • 


* s s 2 V dr o 4 u r o' dr o - u V dr o 4 < u /r o ~ s 3> dd 


+ dt — dt^ . 
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The quantity U may be placed in the form 


where 


U = (t-t Q ) + M ( >p s 4 - 3s^) 


= ip^l 4! + Q'4’^ > /6!+Q!^ tp ^ / 8\ + a^ /10 ! + . . . 

s = ip ^ / 5! + a ip ^ 7 l + a ip ^/9\ + a^ j^^^/11! + .. 


by using the identity 
9s n 


2 —-- = ip s . - 3s 

9a ^4 5 

treated in Section 3.1. 

Thirdly, 

(g d f - f d g) = g(-fdr 0 ^r 0 ~fdr/r- M ds l /rr 0 - d(i s^/ r r Q ) 

-f (s Q dr Q - g dr + r Q ds Q + (T^Sj + d cx^/r 
= (-f g/r Q - f s Q /r) dr Q - (f r Q /r) ds Q -(gM/( r r o) + f^/rJdSj 


- (f s j / r ) d CT 0 - (g s j/ r r Q ) d M 

f r. 


9s, 


= (-f g/r Q - f s Q /r) d r Q - —^ ( <* Sj dip + -— da) 


9s, 


9a 


-(gM /(r r ^)+ f a Q /r) (s Q dip 3—da ) 

-(f Sj/r) d a Q - (g s^r r Q ) d/i 

= (-f g/ r Q ~ f s Q /r) dr 0 + f~ f ^“Sj/r 2 

f r 9 Sq 

- (gM/r r Q + f a Q / r ) s Q /rl r# - [ 

9s i 

+ (gM/r f ^/r) —— ]da - (f Sj/r) d ^ " ^ s 1 / rl () ) d ' 1 
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But the coefficient of r is 


H S 1 r 0 asj 


- f r 0 “Sj/r -(g^ /r r Q + f 

z r ^+ a s, as u \ 

_ / 0 0 0 1 jx _ H 1 _0_ \ 

V r r r rr r ' 


r r r 

0 

<7 


< s o " as i 2) 


and the identity 

2 2 

S — a s 
0 1 


= 1 


is treated in Section 3.1. 


Thus 


(g d f - f d g) = 


( - f « /r o “ £ s o /r) dr o 


3s. 


9s 


-ii- 


[- (r 


0 9a 


0 9a 


3s 


+ M 


9a 


) da 


- s ld r o -s 2 d a Q _ d M + dt - dt Q l 


f r 3 s 

-[—2- —2. + 

r 9a 


g M / r Q + f <r Q 9 s 


3a 


2l da 


-(f s 1 /r) d - ( g Sj/r r Q ) 


dM 


t- fg/r 0 " fs Q /r + 


9 s, 


(r 


0 8a 


t S 1 1 dr 0 


9s, 


+ a 


0 9a 


+ M 


8s 3 

TET ) 


fr 3 s 

0 0 


3a 


g M / r 0 + f ff ( 


3a 


da 


37 



+ ( 


S 2 


-r-» - f 


s„ - 


§ s , 


r r 


) d/i 


— (dt - dt Q ). 


But the coefficient of da is 


T (r o 


3s i 

-— + O’ 

9a 0 


9 s 


r s + cr s 

I ( 0 0 0 1 


9s, 


9a 


+ M 


9a 


M s 1 


-) + 


r r 


M s, 


0 


9a 
9 s 


o’ ) 

0 


1 


9a 


M r 

T [r o 


and the identity 

2 ( s, 


9 s, 


9a 


+ o’. 


9 s, 


9 s. 


0 


9a 
9 s 


1 


9 s, 


+ M 


1 9a 0 9a 

treated in Section 3. 1 below gives 


9a 


+ r (s, 


> = S 1 S 2 + 5 


as. 


8a 


n Hi 

3 r 0 8a 
r 


9 s i 9 s 

+ cr 


0 9a 


2 3S 3 S 1 S 2 +S 3 

+ M —7-- + r 


9a 


] 


8a 9a 2 ^ 


s, s. 


3 to — 

r 8a 


0 




S 1 S 2 d5 l S 0 S 3 - “2 

r -4-=- + r n < —- + -M-) + CT n (- 

2 0 0q/ 2 0 Qqj 


9 s„ 


S 1 S 3 . 
+ —T“) 


s s 9 s o 

2 3 ° 3 1 

+ * -r- + " “sr 

for the coefficient of da . Also, the identities 


9s 


1 


+ s s = s s„ 
9a 0 3 12 


8s 


2 ^ 2 
- + s s = s 

9a 13 2 


treated in Section 3. 1 give 
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t ,r Vz + W % S Z + 'V ^'- 2 IT' ^ 


^ s z s i 

- — + -^~r~ U, (r s, + (T S + n s^) + g -2 

r r 3 201 0 2 4 

r 


L 3 / 2 


- (g- 1) Sj/r + -iy [s 2 (t-t Q ) + g • 2 


= [ —(g— IJSj/r + 


3 Ul /2 


for the coefficient of do- . Thus 

(g d f - f d g) = [ -f g/r Q - f s Q /r + s ^ dr ( 


+ [ -(g-1) Sj/r + -ty U ]( 


0 — — dg . 

— + r • dr - - ) 

2 0 0 r 


fs i gS l 

V s 2 - — > d % + ( Jt T s 3 - TT > dd 

r r 0 


3 (dt " dt 0 ) 


= [-f g/r Q - f s Q /r + s r (g-1)- 2 + Hr “I” U 1 dr 0 

r r r o r r o 


+ ( -(g-IJSj/r + -y-Ul7 0 * dr Q + ( 


M a fs 

—T--- ) d 

3 r 0 

r 


M+(g-D-rr - j- -r + V » 3 - 77 

0 r 0 r 0 


-v < dt - dt „) 

r 
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[-f g/r Q - f s Q /r + ^3 Sj + (g- 1 ) f/r Q + 

(g-D s ! 


+ l - 


3 2 U J dr O* r 0 /r O 

r r o 


3 U1 r 0 ' dr 0 + 


(- 


(g-n fsj 


-) d ff 

J 0 


+ ( s- 


r r Q 


S 3 ) ]d<i " ~3~ (dt - dt 0 ) 

r 


= [ - f s / r r + 


0 0 r r 


+ [- 


(g- 1 ) s x 


0 r 


— 1 v d V ( - ,,r o ,r o 3) - d:r o 

r o r 


f s + (g- 1 )/r 

-Au -i-- '=• 


0 0 


(r * dr + r•dr„) 
' 0 0 0 0 


+ [ - 

or ,finally > 

(g d f - f d g) = 

+ - 


r r 


0 


-f(—=- + 7 

r r L 

0 r 


(g- 1 ) s l 


—-s 3 ) ] dg - (dt - dt 0 ) 

r o r 


1 1 


— + —> *o * dT o - "i" u? 0' dT 0 


i JV f s j + (g-U/r . 


3~ U 1 r 0 dr 0 * 


(r „ • dr + r • • dr ) 
0 0 0 0 


+ [ - 


JL- ( -JL 


- s_) 1 dn - -f- (dt - dt ) . 


Lastly, since 

f i - f g = i 

it follows that 

d f-g + f d g - df*g - f dg = 0 

(- gdf + fdg) = -(g d f - f d g) 
f s, + (f- 1 )/ r 


r • d r„ + f s r • d r„ + f s„ r • d r 
00 200 20 0 


♦ 

+ (g-i) s 2 r 0 *dT 0 + (s 2 /r 0 )dM . 


40 



Then substitution of the final expressions for (g d f - f d g), 

• * * • * 

(-g d f+ f d g), (g d f - f d g), (-g d f + f d g) into the expressions 
above for dr and dr gives the final expressions for the differentials. 
Firstly, 


cfr = £ dr Q + gdr^ 7 


f s + (f - 1) / r i 

1 0 _ , ^ U 

- r * d r - f s^ r * d r — f s_ r * dr 

r _ 0 0 ZOO 200 


0 


- (g-1) s 2 r 0 -dT 0 - (s 2 /r 

r (f-D 8, 


o' <"] 


+ r 


—— V dr o + (f-1) s 2 r Q .d t Q + (f-1) s 2 t Q -d t Q 


♦ * 


+ gS 2 W^O^O - U V d V <U/ V S 3 ! d ^ 

+dt - dt ] 

°J 


or, as the final differential dr, 


dr=fdr+U rr * dr 

0 0 0 


f Sj + (f-l)/r 0 


(f-1) 


r r 0* dr 0 " f S 2 r r 0* dr 0 + T 


s i r V dr o 


t « 


+ (t -" s 2 r V d V « dr o~ u rr o' dr O 


• * 


- l S 2 r V dr 0 - (B -" s 2 r r 0' dr 0 + (f -" S 2 ' V dr 0 * 8 *2 r7 0 d7 ( 


+ [7 (-s 2 /r Q ) + f ( U /r Q - s 3 ) ] d^i + ?dt - ? dt Q 
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Secondly 


d? = f d? + g dr” 4 "r 

0 0 


• 0 1 1 , — u ii _► 

- )r o dr o- -^- Ur o' dr o 

0 r r 0 


(g-U .j ^ 


r 'dr 4 —— U r • dr 
r 0 0 r_ 0 0 


f s 1 4 (g- l)/r 


(r 'dr 4 r *dr ) 
0 0 0 0 


~t <7 


u 


r r 


- « 3 > df - -*5- (dt - dt 
0 r 


+ r 


f Sj 4 (f- l)/r Q 


V dr 0 + f3 2 V dr 0 + fs 2 r 0’ dr 0 


4 (g-i) s 2 ? 0 . d ?o + (s 2 /r o ) d M 


= fdr Q 4U(-ny 
r 


) 'dr - ~ + ~ ? + •yjr'r' . dr“ 

0 0 V r 2 2/ 0 ( 


0 r 


fs 1 4(g-l)/r ^ £ s i 4 (f-l)/r Q ^ 

- ?"? *d? 4 - ~r r” . dT 

r 0 0 r 0 0 


4 f s 2 ^? 0 *d? 0 4 g d? Q ) ? 0 -d? Q 


(g-i) s 2 _ ^ u f s i + (g-l)/i 


r r *dr A — 
0 0 


r r *dr 

0 0 


• • 


4 f s 7r 0 . d r 0 4 (g-l) s 2 r?. d r 4[T(- —) 

0 

• S 2 "? U ”r 

+ ? — - #» “ 3 ( — - s ) 1 dp- „y(dt - dt ) 

r 0 r ° 
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or, as the final differential cfr, 


* t i ^ »; ^ 

dr = f dr^ +U r r^'dr^ 


- i i f s +(g-l)/r_* 

-K~r + ~2 + "2» T? 0' d? 0-7- ' V 10 

0 r r o 

f ®i f (f- 1 )/r , • 

+ --- T v d7 o + f 8 2 7 V d? 0 + g d7 o ■ uT V dr o 

0 

_ ? ? 0 . d ? 0 - li^!i ? r. «* 0 ♦ i s 2 ?? 0 -d ? 0 

r r 


+ (g-l) s 2 ?f 0 .d7 +[ ?(— s t / r r Q + ? (s^r^ + ? ( u / r Q “ s 3 > 1 ^ 


• • •* 

+ 7 dt - r dtp 
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2. 3 THE PARTIAL DERIVATIVES 


The differential relations obtained above may be expressed in 
terms of a set of partial drivatives as follows. Let x, y, z be the 
components of r and x, y, z be the components of "r so that x^, y^, 
z Q are the components of r^ and x^, y^, z q are the components of 
Then the differential relations above may be expressed in the 
following matrix notation. 
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g 0 0 dx Q 

0 g 0 dy Q 

0 0 g dz 

— 1 V- U -J 

f Sj + (g-l)/r 
r 







The coefficient of dt in the above equations is a check on the solution 
of the differential equation. The coefficient of dt Q is the negative of that 
for dt which indicates that an increase in dt^ is equivalent to a decrease 
in dt and vice-versa. In most applications of the differential relations, 

• t • 

t and t are treated as fixed quantities and x, y, z, x, y, z are con¬ 
sidered as functions of x^, y , z^, x^, y^, z^ and sometimes ^ as well. 
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The differential relations are usually expressed in terms of a set 
of partial derivatives as follows. 



8x/9y Q 3 x/0z q dx 0 

3y/3y 0 3y,'9z 0 dy Q 

3 z/ 3y 3z/ 3z dz 

0 0 _ 0 J 


3 x/ 3 x q 3x/ 9y Q 9 x/9z Q dx Q 

+ 3y/3i 0 3y/3y 0 3y/3z Q dy Q 

3 z/ 3x 9 z/ 9 z/ 3z dz 

'■* u 0 0 J L 0- 




3x/3y Q 
3y/3y Q 
8z/ 3y Q 


3x/9z o 
3y/3z Q 
9z/9z^ 



9 *' 9^0 3 */ 3y Q 3 x/ 9 z q "j dx 0 " 

+ 3y/3x Q 3y/3y Q 3y/3i Q dy Q 

_3z/3x 9 z/ 9y 8z/3z di . 

u u 0 J 0 
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The partial derivatives of x, y, z and x, y, z with respect to p 
are merely the coefficients of d^ in the explicit differential expressions 

above. 


8 x/ 8 /lx 
9 y/ 8/1 
8 z/ dfi 



8x/ 8pt 

ay/ dp 

0z/ 8/x 



The partial derivatives of x, y, z and x, y, z with respect to 
x v z and x , y , z are given by collecting the coefficients of 

0' y 0’ 0 00 o 

dX 0* dy 0’ dZ 0 and dX 0* dy 0’ dZ 0 in the explicit differential expressions 
above. A more compact matrix notation is also used to express the 

results as follows. 
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9 x/ 9 x q 

9x/9y 0 

9 x/ 9z 

1 

9y/9^0 

9y/9y Q 

9y/ 9z 

9 z/9 x q 

9z/9y 0 

9 z/ 9z 


g 0 0 x 

0 g 0 -U y 
0 0 g z 


x o y o 
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In the expressions above for the partial derivatives, the parameter 


U = s 2 (t-t Q ) + V- (<P s 4 - 3 s 5 ) 

is a monotomcally increasing and unbounded function of the time t. In 
the case of elliptic motion where a is negative, all other quantities in 
the partial derivatives are periodic with the period 

T = 2 * M y -a 3 

of the elliptic orbit. In general, the partial derivatives of x, y, z , 

• . • • 

y, z with respect to x^, y^, z^, x^, y^, z^ and ^ then have an un¬ 
bounded oscillation in time since the velocity components k, y, i an d 
acceleration components 'x, y, z which are multiplied, by u have an 
oscillation about zero with the period T. However, this effect can be 
nullified for some derivatives if some of the other components x , y , 

Z 0 and V V *0 mult iplying U are identically zero. In particular, the 
coordinate reference frame may be chosen so that only x and y are 

not identically zero, and the unbounded oscillation will then occur in 

only eight of the thirty-six partial derivatives of x, y, z , k, y, i with 
respect to x Q , y Q , z Q , x q , y Q , z Q . 

In the case of elliptic motion, it is shown in Section 3. 1 that 


^hich, with 


A E/2! - (I-cos AE) 
——— - 

\fin a 

3 

A E/3! - ( A E - s i n AE) 
~~~ ~ 5 “ 

v/m /a 


<p = 


AE 


S 2 = 


\/m / a 

1 - cos A E 


M /a 


f \ „ sin A E 

{t_t - ) - r o + % 


g = 


v/fTTaT 

sin A E 




+ Q 


I- cos A E 
M / a 

1 - COS A E 


A E - sin A E 

\Z~V~Ta 3 


M /a 


50 




1 - cos AE 


sin A E 


from above, gives 
AE 


U = 3 \x 


sin A E 


/ M 


n /a 


(g - M 




Of course, the partial derivatives may be computed for elliptic motion 
by using AE rather than >p to solve for the coordinates and determine 
all the quantities in the partial derivatives. Then the only secular term 
in all the derivatives is that due to A E in the expression above for U . 
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2.4 THE INVERSE PARTIAL DERIVATIVES 


The solution of the two-body problem is still of course valid if 
t is treated as the reference time and *r and ~? as initial conditions. 
Then r^ and r^ are obtained as solutions of the differential equation 
at the time t^. For this inverse solution, 


r 



a = ~r • ■? 
a = *?•"? — 2 ^ /r 

where a is constant in time and therefore equal to the value above 
computed from if and 'f . As in the solution above, 1JT and its 
transcendental functions. 

8=1 + a if/Zl + a 2 4 /4! + a 3 f 6 /6! + ... 
s = "$ 4 ar ^ 3 /3! + Qf 2 ^ 5 / 5 ! 4a 3 ^ ? /7 ! 4 . .. 

S 2 =^ / , / 2 ! 4 a; i ]) / 4 ! 4 a $ /6!4{* 2 $^/8!4.... 

s 3 ="? 2 / 3 4 a j/T 5 /5!4 a 2 ^ 7 /7!4Q! 3 ^r 9 /9!4... 


satisfy the equation 

^ = t 4 r Sj 4 ffl 2 4 n l y 

The bar is used to distinguish parameters of the inverse solution 
from those of the solution above. Also 

r 0 = r V-l + “ 7 2 

and 


f = 1 - fJ. / r 

- “ V r o r 


g= , 0 - , -*= 3 
g = 1 -H = 2 /r 0 
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which give 


r = f r + g r 

0 5 

= Tr + g r 

0 6 

as the inverse solutions of the differential equations. 


However, 


* /•*. 
J r* 


from the definition of ip above where r* is the value of the radius at 
the times t* between t^ and t. But 


t o r 1 f 

/ _dt* = - / dt * 

t J r* t J r* 


so that ip is merely the negative of ip . Thus 


= s 

0 0 


s = — s 
1 1 


S 2 = S 2 


= - s^ 
3 3 


so that ip satisfies the equation 


Also 


V *- ' S 1 + as z - " V 


r 0 = rs 0- <r *l + “ S 2 


f = 1 - M s /r 


g = t 0 - t + „ Sj = - (t-t 0 - m s 3 > 


f = A* Sj/ r r 0 


g = 1 - u s,/r 


2 0 
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This gives the formulas 

r Q = gf - gt 

r = -f r + f *r 

0 

which were used above to obtain the partial derivatives, 
substitution of 


r = f r + g r 

# 0 6 .0 

r 88 f r 0 + g r o 

into the first equation gives 
= (f g - £ g) r n 


In addition, 


so that 

f § - i g = i 

which was also used to obtain the partial derivatives. Another property 
of the inverse solution is that 

U = s 2 (t Q -t) + M ( J 3*4 - 3r 5 ) 

where 



which gives 

U = - U - 

The relations above may be used to express the partial derivatives 
in the inverse differential relations 
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pV 3x 

% = 8y 0 /9x 

-*oJ K /8x 


9x / 9y 3X /az dx 3x /ai 


9y 0 /3y dy Q /dz dy + ay Q /ai 
te Q / dy az fl /9z _dz j _^ Q / 8X 


ax Q /9y 

3y Q /ay 
a z Q /ay 




dx o] pV 8x 
% = 9 y 0 /9x 

_ d *oJ l a V ax 


ax Q /ay dk Q /dz r- 


dz Q /dy az Q /a zj [_dz 


3x^/ 8x 


dy Q /dy dy Q /dz dy + dy Q /dk 


az Q /ai 


ax Q /ay 

dy Q /dy 

dz Q /dy 


ax Q /az 

ay 0 /9z 

az Q /az 
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A similar method of expressing the partials of x Q , y , z^ and 
x , y Q , z Q with respect to x, y, z and x, y, z by first using the 
functions of W and then converting to functions of i/ gives 


9 x /9 x 9 x Q /9y 
9y Q /9x 9y Q /3y 
9 z^/ 9x 9z q/ 9y 


9x q / 9 z 

9y Q / 9 z 

9z q / 9 z 


g 0 0 


o g o -U y [ 


0 0 g 


x y z 


r o-f s + (g-l)/r 
x^ x^ 1 

0 0 - - 

r 

+ y 0 *0 —<*- 1) s i 


x y z 


(g-l)s 2 X y 


9x q / 9 x 9 Xq/ 9y 
9y Q / 9x 9y Q /9y 

9z / 9x 8z / 9 y 

- 0 0 


9x q /9 z 
9y Q /9z 
9 Zq / 9 z 


g 0 0 x Q 

o g o + u y Q C* y 7 ‘. 

° ° g z n 


x 0 X 0 r £ S 2 (f - 1)s 2 X y Z 


y 0 y 0 L“ ( *“ 1)9 2 gs 2jl x y z 


_ z o Z ( 


9 x Q / 9x 


9x / ay 9x / a: 


9y Q /9x 9 y Q / 9 y 9 y Q / 9z 
9z^/ ax 9z rt /9y 9z^/ 9z 


z l = - 


f 0 0 

0 f 0 -U 


y o L x 


0 0 f 


-f(— + Y + 

rr o * 

f Sj + (g-l)/r 


f s 1 -t (f-l)/r Q 


y z 
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0 


~3 x Q / & 

9x Q /3y 

3^/3^ 


f 

0 

0 


*0 

3y 0 / 9 i 

3y Q / 3/ 

3y Q / 3z 

= 

0 

f 

0 

+ u 

^0 

3 z Q / 9x 

3z Q /3y 

Bz^ / Bz : 


0 

0 

f 


• • 

~ z 0- 


+ 



Yr 



f Sj + (f- 1)/r Q (£-1) s 1 


-is 2 (f-D s 2 J 



These inverse partial derivatives are often utilized for computation. 
However, they satisfy identically the well-known formulas 


B x / Bx 

0 

3x Q /3y 

Bxq/Bz 

3y 0 /3x 

3y Q /3y 

3y Q /3z 

a Zq/ Bx 

3z Q /3y 

8Zq/ 3z 

B x^/ Bx 

3 x Q /3y 

3Xq/3z 

3 Y{) / 3x 

3y Q /3y 

ay Q /az 

3z / ax 

V- (J 

3z 0 /9y 

a Zq/ az 

ax Q / ax 

3x q / 9y 

3Xq/ az 

3y Q / 3x 

3y Q /9y 

ay Q /8z 

3 z^/ 3x 

3z Q /ay 

3Zq/ az 


Bx/ Bx^ 

3x/3y o 

3x/3 z q 

ay/3 Xq 

3y/3y Q 

3y/3 Zq 

_3z/ 3Xq 

3^/3y o 

3z/ 3 Zq 

~3x/ ^Cq 

3 x/3y Q 

3x/ 3 Zq 

3y/^X q 

ay/3y o 

3y/ 3 ^q 

3z / 3 Xq 

az/ 3y o 

Bz/ Bz 

0 

~ax/ ac Q 

3 x/3y Q 

Bx/ Bz q 

3y/ 3 Xq 

3y/3y 0 

3y/3z 0 

az/8 Xq 

3z/3y 0 

3 z/8 Zq 


9x q / 3 k Bx q / By 8x Q /az 

By /Bx By / aJ" By /B z 

7 0 7 0 7 0 

B z / Bx Bz / By Bz / Bz 

0 0 0 -J 


ax/ax Q 

3x/3y Q 

3x/ 3 Zq 

T 

3 y / 3x q 

3 y / 3y 0 

3 y / &z 0 


3z/ 3x 

L. 0 

3z/3y Q 

B z/ Bz 

0 J 



where the T indicates matrix transposition. Thus the inverse partials of 
Xq, y^> z^ and x^, y^, z^ with respect to x, y, z and x, y, z may be 

t # # 

obtained directly from the partials of x, y, z and x, y, z with respect 


to x 


o’ V 


z and x . y , z, 
0 0 y 0 ( 
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3.SOME USEFUL FORMULAS FOR DERIVATION AND COMPUTATION 


3. 1 THE PROOF OF SERIES IDENTITIES 

Several identities were used in the derivation above for the partial 
derivative s. 

As an example, consider 
2 2 i 

5 o - " 5 1 '• 

This may be proven by collecting coefficients of the same powers of $ 
in the products of the series. 

(1 + a 4 2 /2! + a 2 //4! +_) • (1 + a ip Z /2 \ + « 2 ip 4 / 4! +-) 

-a (4 + a i/- 3 /3! + a 2 4^/5 !+ ...) ( j/j + ck ^ 3 /3! + « 2 ^ 5 /5! + ....) 

= 1 + a 4 Z { 1/2 ! + 1/2 !) + a 2 4 4 ( 1/4 ! + l/2 ! 2 l) + 1/4!) + .... 

- a 4 - a Z i/> 4 (l/3! + 1/3!) 

= 1 + a^ 2 (l-l) + a' 2 i^(l+6+l-4-4)/4! + .... 

= 1 


The product is valid since the series s . s . s_, s , s , s_ above are 
r 0 1 Z 3 4 5 

always absolutely convergent. However, an alternate proof is to first ob¬ 
serve that the identity is true for zero a where each series reduces 
to its leading term. If a is negative, then 

a = - [i I a 


S q = COS AE 

Sj = sin AE 
vV /a 

From above where 

AE = sj /u/ a 4 

so that 

2 Z Z _ u 

s_ — as, = cos AE 4 r 
0 1 a 

= 1. 


sin 


s/VF< 


AE 2 


= cos 


AE 4 sin AE 
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Similarly, if a is positive, then 


s Q = cosh AF 


s j = sinh 1 AF 


from above where 

AF -\/a i 

so that 


2 Z Z 

Sq - a = cosh AF - 


= 1. 


sinh. AF .2 Z Z 

-) = cosh AF - sinh AF 

vAT 


Thus the identity is true for all a and all $ . 

The other series identities used in the partial derivatives involve 


3s 3s, 

0 1 


3s, 


3s. 


They may also be proven by collecting 


da * da 9 da ’ da 
coefficients of the same powers of ip in the products of the series. 

However, the relations 

3 s 

0 = *s 
c 0a 1 


0S 


da 

3s^ 


i- = * S 2 - S 3 


-r- = Ip s Zs 

da 3 4 


3s, 


da 


= ip s 


3s, 


may be proven by a term-by-term comparison Of the two sides of each 
equation. If these formulas are substituted in the identities, each identity 
is valid for zero a . But 
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s 0= l+a» 2 

s = a s 3 

s 2 = 4 > 2 / 2! + a S 

^ 3 

s = ip / 3! 4 as 
3 D 

similarly, so that, for non-zero a , 


s o “ 1 


1 - s 


0 


2 a - a 

cos h AF - 1 1 — cos AE 


a 


s 3 


a 


/i/a 
- S 1 


- a 


sin h AF/ \f~6T - AFJV aT AE/ \/ jj. / aT - sin AE/ \/ /x /a 


a 


/i/a 


sinh AF — AF 


S 4 = 2 


AE — sin AE 

/oT 3 , v/fiTa 

s 0 - * /2! * /2! - s 0 


a 


a 


(cos h AF -1)/ g - (AF/v¥) 2 /Z !_ ( AE/!-(1 - cos AE)/( £/a) 
a /i/a 


(cos h AF -1) -A F/2! 


a 


A E/2! - ( 1 - cos AE) 

.2 


( « /a) 


s - ip / 3! 

S 5 '-A- 

a 


<T/3! - s. 


a 


(sin h AF - AF)/ /?"-( AF//lT) 3 /3 !_ (A E^TT^) 3 / 3-( A E-sin AE)/ft /g 3 
a /x/a 

(sinh AF - AF) -A 3 F/3! AE/3! -( AE - sin AE) 


/a"' 


sfVTi 
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Use of these relations proves the identities for both negative 


a where AE is real or positive 0! where AF is real. For example, 


3s 


S 2 S 3 = 2 (S 2 


1 


3a 

3s. 


1 3a 


= s z ( Z- 


“—) - s (2 

3a 1 


3s„ 


3a 


s^ ( s z - s 3 ) - s l ( 0 s 3 - 2 s 4 ) 


is true when a is zero since 

{ip 2 /2!) ( ^ 3 /3 !) = ( 0 2 /2!)( >p 4> 2 / 2! - <p 3 / 3!)-( 4> )( <l> * 3 /3 ! - 2^ 4 /4!) 

i/) 5 / 12 = i/) 5 (1/4 - 1/12 - 1/3 + 1/12) 

= ip / 12 . 

For negative a , 


1 - cos AE AE - sin AE 

~ _ ^ 

v/ /i/a / M /a 


1-cos AE 
M / a 



j.— c o s £E 
M/a 


AE - sin AE \ 
v/m /a 3 / 


sin AE j AE AE - sin AE 
v/m /a VM/a \/m /a 3 


- 2 


A E/2! — (1 -cos AE) 

( M / a ) 2 


(1-cos AE) ( AE - sin AE) 


(1-cos AE) (AE -AE cos AE - AE + sin AE) 

2 

- AE sin AE( AE - sin AE) + 2 sin AE A E/2! 

-2 sin A E (1 - cos AE) 

(1 - cos AE) (sin AE - AE cos AE) 

2 

+ AE (1 - cos AE) - 2 sin AE (1 - cos A E) 

(1 - cos AE) (sin A E - AE cos AE + AE + A E cosAE 

- 2 sin aE) 

(1 - cos AE) (AE - sin AE). 
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A similar result follows by using AF for positive a, but the use of 
AE is actually sufficient since the trigonometric formulas used are 
valid for imaginary arguments. 
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3. 2 THE USE OF DIMENSIONLESS SERIES 


For many purposes it is convenient to define the dimensionless 
argument 


X = a ip 

and its transcendental functions 

c = 1 + A/2! + A 2 /4! + A 3 /6! +- 


c = 1/1! + A/3! + 

c = 1/2! + A/4! + 

A Z /5! + 

A 3 /7 ! +- 

A 2 /6 ! + 

A 3 /8! +- 

c. 

c_ = 1/3! + A/5 ! + 

A 2 / 7 ! + 

A 3 /9 ! + .... 

3 

c = 1/4! + A/6! + 

A / 8 ! + 

A 3 /10 ! + .... 

4 

2 , 

3 , 

c = 1/5! + A/7! + 

5 

A / 9 ! + 

A / 11! +- 

which are related to s^, 

V V 

S 3’ S 4’ S 5 by 



or, conversely, if A ^ 0, that 
C 5 = (C 3 " 1/3!)/ * 

c 4 = < c 2 ~ 1/2! )/ A 

c 3 = (c 1 - 1)/A 
C 2 = <<=0 - !)/*■ 
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In terms of trigonometric or hyperbolic functions, for A ^ 0, 


c Q = 1 + X/2! + A /4! + A /6! + 


2 4 6 

= 1 +/T /2! + /T / 4 ! +/T / 6 ! +_ 

= cos h /*= cos h(i v/^A) 

= cos \/ - A 


c x = 1 + A/3! + A / 5! + A / 7 ! +_ 

3 3 

= ( /~A+ /T/3!+/a/5!+_) //a” 

= sin h/T” //T = sin h (i /7 a ) / (i /77) = i sin /^T/(i /7 a~ ) 
= sin /7T / /-A . 

Let Cq {4A} and c^{4A}be, respectively, the values of c^ and c^ for the 
argument (4A) rather than A . Then, forA>0 and A < 0, respectively, 
c q {4A} = cos h ^ 4A = cos y/-4~A 


= cos h (Zv^ A ) = cos (Z v/ -A ) 

= Z cos h^\/ A -1 = Z cos^ / - A -1 
c^ {4A} = sin h V^4 A / \/ 4A = sin -4A / l/ -4 A 

= sin h (Z /a ) / (Z/T ) = sin (Z \J - \ ) / (Z \J - A ) 

= Z sin cos h /r / (2/7) = 2 sin \J - A cos \/ - A / (Z ^/- A ) 

= cos h/7: (sin h/7//7 ) = cos y— A * (sin /77/ /7 a) 
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so that 


c 0 { 4A } = 2 c 0 2 { A > - 1 
c, {4X}= c Q { X t ' c, {X} . 

These two relations are also true for zero X , and are thus valid 

for all values of X . It is worth noting that the argument 

. ,2 

X = a ip . 

and a are simultaneously positive, zero, or negative except in 
the trivial case where ip 3-nd thus X are both zero. 
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3. 3 FORMULAS FOR THE SOLUTION OF KEPLER'S EQUATION 
It was shown above that the general equation 

' = V r 0 S 1 + 0 S 2 + MS 3 

relating t and * is a generalization of Kepler's equation for elliptic 
motion. Stumpff (1962) calls it the main equation but prefers an 
equivalent of the alternate form 

1 = V r 0 * - °0 s 2 + T o V 

The first, second, and third derivatives of this alternate form 
with respect to </> give the formulas 

r = r + <t s 4 y s 

0 0 1 y 0 2 


for the radius, 


s„ + y s 

0 0 0 1 


for the scalar product of the position and velocity vectors, and 


r * v o s o + “ °b s l 


for the parameter y which equals /i for parabolic motion, i. e. for 
“ = 0* When CT 0 and y 0 are bot h zero, r is always equal to r 
whereas a and y are both always zero. Then the motion is circular 
and is a linear function of t. 

In computing coordinates with or without partial derivatives for 
a particular time t, it is necessary to solve for that value of ip 
which satisfies the main equation above. With a slight change in 
terminology, let 

u-y = y,t yyc 
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be the time interval computed directly from any value of ip , 
and let T be the desired time interval. Then the change in 
ip required to change (t - t^) to t is given by 

T= (t - t ) + d ^ [ V Alp 

° d>P 

to first order accuracy. Thus 
T- ( t - t Q ) = r Alp 

where r is the radius for the particular value of $ . The value 

ip* = ip + A ip 

is then a better approximation for the desired value of $ if the 
first order approximation is accurate. Explicitly, if 


At 


r 0 *1 + S 2 + “ S 3 


T 


r - r o s o* Vi + “ 8 2 


are first computed from then 

ip = ip - At / r . 

This is Newton’s method for the solution of the main equation. 

* 

Using ip as a new value for ip then gives a new approximation, 

etc. However, a method of stopping the iterations is necessary. Also, 

Newton’s method is not always convergent and should be backed up by 

* 

alternate methods for obtaining a better approximation ip from ip . 
The iterations may be controlled by always bounding the solution by 
the maximum known value ip and minimum known value ^ + for which 
the respective residuals At_ and Ar + are, respectively, negative and 
positive. Then any approximation is accepted at a new value of ip 
only if 

ip_ < $*< ip + 

* 

i. e. if 4> is between ip_ and ip + but not equal to ip_ or ip + . 

* 

When the value $ computed by Newton’s method fails to lie 
between and , alternate methods may be tried to obtain a 
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satisfactory value of # . One such formula is 
/ = ip (1 - At/t) 

which gives smaller corrections to for smaller values of |At| . 
Another is to define tp* as the value satisfying the interpolation 


formula 

= 



At_ At + 

Explicit solution of this equation gives 



At 

At + — At 



The use of these formulas is treated in the description of the 
Fortran subroutine in Section 4.1. 
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3.4 FORMULAS FOR SERIES COMPUTATION 


Each iteration in the solution of the general form of Kepler's 
equation requires the evaluation of s^, s^, s^, and from a given 

value of ip . The evaluation of and from the final value of ip 

is also required for the determination of partial derivatives. However 

it is convenient to compute the dimensionless series c^, c^, c^i c^, 

c , c of the argument 

4 5 

2 

X = a ip 

rather than s rt , s . s_ , s . s . s directly. Then 
0 12 3 4 5 


s o " c o 

s i = * c i 


s 2 = * C 2 
,3 

s 3 = * °3 


give the values required for each iteration of Kepler's equation. The 

final value of c and c may be used to obtain 
4 5 


U = s 2 (t-t Q ) 4 M (c 4 - 3c 5 ) ip 


so that s and s r are never required explicitly. 

4 5 

The series c , c . c , c . c . and c c must be computed from 
0 1 2 3 4 d 

arbitrarily large values of M , i. e. for arbitrarily large values 
of |o?|and \ip\ , in order to determine coordinates and partial de¬ 
rivatives for arbitrarily large values of |t|, i.e. for any value of t 
This may be done by repeatedly dividing X by 4 when | x\ > 1 until 
the m division reduces the argument to one or less in magnitude. 

Then c and c, are computed as described below with the reduced 
0 1 r 

argument. Then the formulas 
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c o M = 2c 0 2 (4 -1 

c i = C 0 W ' c jW 


are applied repeatedly m times to obain their true values for the 
original value of X . Then the formulas 



c 3 = (c 2 - 1)/A 

c 4 = (c 2 - 1/2) IX 

Cg = ( c 5 - 1/6) /A 

above give the true values of the four remaining .series. This approach 
seems to be quite satisfactory, but an alternate method is to determine 


Cq = cos V-A 

0 ^ = sin v r - A / /“A 
directly when A< - 1 or 
Cq = cosh ^A 

Cj = sinh aTX / v^A 

directly when A > 1 . Then the same four formulas above give c , c , 

2 3 

and c . 
b 

The series c Q , c^ c^, c^, c^, c^ may all be determined from a 
given value of A by first using the nesting formulas 



— 

A 

X 

\ A ~ 

(2n+3)(2n+2) 

(2n+l) (2n) 

/ 9*8 

X 

X 

. . . A A ~ 

(2n+2)(2n+1) 

(2n)(2n-1) 

J 8-7 

J 




to determine c^_ and c . The number of terms n required for the 
computation is considered below. Eight terms, i. e. seven leading l’s, 
are sufficient for accurate floating-point computation to sixteen decimal 
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digits when |x|< l which is always true since if | x| > 1 initially 
it is either reduced as described above or the alternate method is 
utilized. Then the resulting values of c^ and c^ give 


c^ = 1 / 6 + Xc^ 
c~ = 1/2 + A c . 


C 1 = 


C 0 = 


1 + X c„ 


1 + X c_ 


for the remaining four series. If the argument has been reduced, c Q 
and Cj are then used as c 
c , C-., c_, c _. c., and c. 

A ' 1 J ' v A !■ 


and c^ are then used as described above to obtain the true values of 


'1 


'5. 
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3. 5 THE ACCURACY OF THE SERIES COMPUTATION 

The limiting cases in the nesting computation above occur when 
| = 1. If -1, then 

Cq = cos V -A = cos (1) 

= .54030 

C 2 = (c Q - 1) / X = (. 54030 -1)/ (-1) 

= .45970 

C 4 = (c 2 - 1/2) / A =(. 45970 - .50000) / (-1) 

= .04030 
and, similarly, 

C 1 = sin vUT /\/ — A = sin (I)/(1) 

= .84147 

c 3 = ( Cj - 1) /A = (. 84147 - 1) / (-1) 

= .15853 

c 5 = (c 3 - 1/6) /A = (. 15853 - . 16666) / (-1) 

= .00813 

Also, if A = 1, then 

c Q = cosh VA = cosh (1) 

= 1.54308 

C 2 = {c 0 " 1)/A = f 1 - 54308 - 1 )/ (!) 

= .54308 

c 4 = ( c 2 - V2)/A = (.54308 - . 50000)/(l) 

= .04308 
and } similarly, 

c 2 = sinh = sinh (1)/(1) 

= 1. 17520 

C 3 = (c i " X)/ A = O. 17520 - 1)/(1) 

= . 17520 

C 5 = (c 3 ' 1/6)/ A = (• 17520 - . 16666)/(1) 

= .00754 
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Suppose that eight terms are chosen to compute c 4 and c 5 for 
(a(= 1. Then if A = -1, the truncation error Ac^ introduced in 

computing c is 

4 8 9 10 

A c = A / 20 ! + A / 2Z ! + A / 24! + .... 

4 


= + 1/20! - 1/22! + 1/24! + - 


so that 


Ac < 1/20! = . 41103 x 10 

4 


-18 


from the properties of alternating series. 
Thus 


Ac -18 

^"4 < .41103 x 10 


04030 


= 10. 19 x 10 


18 


Ac, 


< .1019 x 10 


16 


Similarly, the truncation error Ac,, introduced in computing c^ is 

8 9 10 

Ac — A /2 1 ! + A /23 ! + A /25 ! 4 • • • • 

5 


1/21! - 1/23! + 1/25! - ... 


so that 


A c < 1/21 ! = . 01957 x 10 

5 


-18 


and 


Ac 5 < ^01957 x 10" 18 


.00813 


= 2.407 x 10 


-18 


Ac 


— | < . 2407 x 10" 17 
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If A = + 1 with eight terms, the truncation error Ac introduced 

4 

in computing c^ is 

8 9 10 

Ac 4 = A /20 ! + A/22!+ A / 24!+.... 

= 1/20! + 1/22 ! + 1/24! + . . . . 

< 1/20! + l/(20 2 • 20!) + l/(20 4 • 20!) + _ 

= (1 + 1 / 20 2 + 1 / 20 4 +_)/ 20 ! 

= [1/(1 - 1/20)) / 20 ! = (20/19)/20! = (1 / 19 !)/ 19 

= .04326 x 10" 17 

so that 



.04326 x 10 
.04308 


-17 


= 1. 004 x 10 


-17 


I -1 < . 1004 x 10 ^ 

C 4 

Similarly, the truncation error Ac introduced in computing c is 

^ 5 

8 9 10 

Ac 5 = A / 21! + X/23I+ A /25!+.... 

= 1/21! + 1/23! + 1/25! + _ 

< 1/21 ! + 1 /( 21 2 - 21 !) + 1 /( 21 4 . 21 !) +_ 

= (1 + 1 / 21 2 + 1 / 21 4 + . .. . )/21 ! 

= [ 1 / ( 1 - 1 / 21 )] / 21 ! = ( 21 / 20 )/ 21 ! = ( 1/20 ! )/20 
= .02055 x 10' 18 
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so that 


A c — 18 

5 < . 02055 x 10 


.00754 


= 2. 72 x 10 


-18 


< .272 x 10 


Thus eight terms are sufficient for sixteen significant figures 
in the nesting computation above for c^ and c 4 with |a|< 1. A 
similar derivation will determine a sufficient number of terms if a 
different number of significant figures is required. Ordinarily this 
is determined by the particular electronic computer utilized, e. g. 
the IBM 7094 has 54 binary bits in a double-precision fraction which 
gives slightly more than 16 significant decimal digits for floating¬ 
point computation. 
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4. THE DOUBLE-PRECISION FORTRAN 4 IBM 7094 SUBROUTINE 


4. 1 DESCRIPTION OF PROGRAM INPUTS AND OUTPUTS 

The Fortran 4 IBM 7094 double-precision subroutine TWOBDY 
has the following nine inputs. The six initial position and velocity 
components at the reference time t are input as a state vector 



The seventh input is the time interval 
T = t — t 

0 

between the reference time t and the time t at which a solution is 
desired. The eighth input is the constant fi in the differential equation 
where 

M = G (m^ + m 2 ) 

if relative coordinates between two gravitationally attracting masses 
m i and m^ are being obtained. The ninth input is an approximate 
value for the final solution i p of Kepler's equation where 

4 > = o 

is used if no approximation is available. 

The input value of >p is changed by the program to the actual 
solution of Kepler's equation for the input time interval r . Thus 
the ninth input to the program is also an output quantity. The number 
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of iterations required to solve the general form of Kepler's 
equation can be considerably reduced if a good input approximation 
for ip is available. In particular, solutions may be required for a 
whole series of values of r which differ by given increments. Then 
the use of the previous value of ip as an input approximation for ip 
at each new incremented value of T is very efficient. It is usually 
better to use an input approximation for ip even if it is only known 
very roughly. If not, it is more efficient to use an input value of 
zero although the program will work with any input value for the 
parameter ip . 

The program outputs the six position and velocity components 
in the form of a state vector 



for the time t. The 36 partial derivatives of x, y, z, x, y, z with 
respect to x , y , z Q , x Q , y^, are output in a 6 x 6 matrix 

ax/ 9 x Q 3 x/ 3y Q 9 x/ 8 z Q Sk/dk Q 3x/3y Q dx/dz Q 

sy/ 3* 0 dy/9y/ 8z 0 9 y/ 8 * 0 d y/ d y Q a y/ 8z 0 

3z/ 8 Xq dz/dy-Q 3z/ 3z^ dz/dx^ &z/ 3y^ 3z/9z^ 

3 x/ 3 x^ 3x/3y^ 3x/3z^ 3x/3x^ 3x/ 3y^ 9x/9z^ 

3y/3x Q 3y/3y Q 3y/ 3z q dyl 3x Q 9 y/‘V q 9 y/ 

3z/3Xq 3z/3y^ 3z/3z^ 3z/3x^ 3z/3y^ 3z/3z^ 
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Similarly, the inverse partial derivatives are output in a 6x6 matrix 


8 x^/ 8x 

9x q / 9y 

9x^/9 z 

9x^/9x 

3x q / 9y 

9 x Q /9z 

9 y Q /9x 

9y Q /9y 

9y Q /9z 

9y Q /9x 

9y Q /9y 

9y Q /9z 

8z^/8x 

9 Z ()/ 9y 

0z^/8 z 

9z^/9x 

9z q / 9y 

9z^/9z 

9x^/8x 

9x q / 9y 

8Xq/ 8 z 

9x^/9x 

9x q / 9y 

9Xq/9z 

9y Q /9x 

9y o /9y 

9y Q /9z 

%/& 

9y Q /9y 

9^ q /9z 

9z q /9x 

9z Q /9y 

0Zq/8z 

cfc^/ 9x 

9z Q /9y 

9^ q /9z _ 


The two 6x6 matrices above are also matrix inverses. That is, their 
product should give a 6x6 identity matrix. 

The six partial derivatives of x, y, z, x, y, z with respect to 
^ are output as a vector 
9 x/ 9fx 
9 y/ 9p 

9 z/ 9p 
9 x/ 9p 
9 y/ 9p 
_9 z / 9p_ 

Similarly, the six inverse partial derivatives ofx.y.z.x'.v.z 

0 *0 0 0 y 0 C 

with respect to M are output as a vector 
9 x Q / 9ju 
9 y Q /9p 
9 Zq / 9p 
9x Q /9p 
9 y 0 / 9p 

_9 z Q /3MJ 
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Additional outputs are the acceleration vector 



«• 

y 


UJ 

at time t and the acceleration vector 



♦# 

z 


L 0-» 


at time t , as well as the 
r = V x + y 
at time t and the radius 


radius 


2 2 
+ z 


r 

0 



+ z 


r 

o 


at time t^. 

The Fortran symbols for the inputs and outputs above are 
listed in comment statements at the beginning of the subroutine. 

The program may be used as a "black box" by those not interested 
in the formulation or computation. However, a knowledge of the 
operation of the subroutine is very important if modifications are 
desired for particular applications. For example, it may be de¬ 
sired to eliminate the partial derivatives or at least make their 
computation optional since they are not required in many cases. 
However, such modifications will be worthwhile only where the 
particular application is encountered repeatedly. Ordinarily it is 
better to use the program as it stands without modification although 
this may involve some redundant computation. In any case, the 
operation of the program is described below for those who are in¬ 
terested. 
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4. 2 THE INITIAL COMPUTATIONS 


The following initial computations are performed only once 
before starting the series computation below that is repeated for 
each iteration in the solution of Kepler's equation. The following 
technique for computing the radius r Q at time t is used to avoid 
utilizing a square root subroutine. First 

b = maximum (| x Q |, | y Q |, |z |) 

c = (x Q /b) 2 + (y 0 /b) 2 + (z Q /b) 2 

are computed and the initial value of d is set equal to two. Then a 
better approximation d# for the square root of c is obtained from 

d* = (d + c/d)/2 

and d* is used as a new value of d to iterate the formula. The 

iterations are continued until d is no longer greater than d* which 

is then accepted as the square root of c. Then 

r = b • d* 

0 

gives the radius r Q at the reference time. The quantities 


i + v V + z z 

0 0 0 y 0 y 0 0 0 

. 2.2 .2 

x o +y o + z o - 2M/r o 

are then computed. It is also necessary to initialized a counter m 
to zero in order to compute the series described below for large 
arguments. 

In the trivial case where t is zero, tp is set equal to zero and 
the series computation below is performed. The known value zero for 
tp satisfies Kepler's equation identically so that all the outputs will 
then be computed by the program. If r is negative, it is known that 
ip must also be negative. Thus the right bound >p for ip is set equal 
to zero and its residual At is set equal 
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to the known value - T , whereas the left bound >P_ and its residual 
At are both set to - °° , i. e. the negative of the largest number 
available for computation. If T is positive, it is known that >P must 
be positive. Thus the left bound j p _ for tp is set equal to zero and its 
residual A T is set equal to the known value -T , whereas the right 
bound i p and its residual Ar + are both set to + °°, i. e. the largest 

number available for computation. 

If the input approximation for 4> then lies between 4>_ and ip + , 
it is used in the series computation below. If not, Newton's method 
using an initial <p of zero gives the approximation 

*= T /r o 

which is then tried. If this value is not numerically between ip _ and 

♦ V 

4 > = T 

is used to start the series computation. Thus an initial approximation 
for i l> is always obtained which is between ip _ and i/» + but not equal 
to either. 

4. 3 THE SERIES COMPUTATION 

The computation of the series is begun by computing the argument 
2 

X = a ip 

of the simensionless series. If | x| > 1> the value is saved as 

but A itself is repeatedly divided by 4 until < 1, and the number 

m of required divisions is also counted. Then 3c & and c ^ are computed 


by the nesting formulas 
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which in turn give 


c 3 = [1/2 + A (3c 5 ) ] / 3 

c_ = 1/2 + Ac 
l 4 


where each series is correct to sixteen significant figures. 

If the original value of A has not been reduced, the series s 

3’ 

s 2’ S 1 are com Puted as described below. If the original A has been 
reduced, i. e. if m> 0, the formulas 

* 2 

c o = 2 c o - 1 

❖ 


are applied m times with c and c as new values for c and c . 

5jC U ft t 0 1 

The final values of c q and are accepted as the True values of c 

and c for the original value of A, i. e. for A . Also 

P 


C 2 = {C 0 - 1}/ % 
c = (c - 1)/ A 

Si P 

C 4 = (C 2 - 1/2 >/ % 
3C 5 = (3c 3 " 1/2)7 X 


are computed as the true value of the dimensionless series. 

The dimensionless series c , c , c give the series s , s , s 

3 2 1 3 2 1 

by the formulas 



S 1 = * C 1 
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which are used to compute the residual in Kepler's equation below. 

The value of c Q is equal to s Q and will be used in Newton's method 
for solving Kepler's equation. If this particular series computation 
should be the last of the iterations for solving Kepler's equation, c 4 
and 3c 5 will be used below to compute partial derivatives. The series 
computation described here seems to be quite satisfactory for ob¬ 
taining solutions for two-body motion. The absolute error in the so¬ 
lutions will increase as T and thus X increase, but this is an 
inherent error since the absolute error in i itself increases as it 
becomes larger if it is represented as a floating point number. 

4.4 THE SOLUTION OF KEPLER'S EQUATION 

The residual At for solving Kepler's equation is found by first 
computing 

8 = r 0 S 1 + % S 2 
for the current value of 4> . Then 

At = g + n S 3 -t 
gives the residual. The radius 

r = r 0 C 0 + Vl + MS 2 

for the current value of ip is also determined. If A T is zero, the 
coordinates are then computed as described below since Kepler's 
equation has been solved. If At is negative, it replaces the old 
value of At and i p replaces the left bound ip_. If A t is positive, 
it replaces the old value of At + and ip replaces the right bound ip + . 
Then Newton's method 

i p = ip - A T /r 

is used in an attempt to obtain a better approximation for tp • 

If ip* > ip and ip* < 4> + , 4>* is accepted as a new value of 4> 
and the series computation above is repeated. Otherwise, the following 
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methods are used succesively in an attempt to obtain an approximation 


between if> and ip . First 

4> = *l>_ (1 - 4 At / T ) 

is tried if |at - |<| At+| but 

4> = f (1 - 4A t /t) 

+ + 

is tried if I At- I >1 At + I 


The factor 4 helps to give significant 


corrections to ip, in floating-point computation when |at - | or |at + | 

is small. If | At-I = |at + | or / fails to lie between 4> and 0 , 

the value 


ip* = 2 >P 


is tried if t > 0 , but 
ip* = 2ip + 

is tried if t<0 . This insures that the solution will always be isolated 
between finite bounds. If this method fails, the interpolation formula 
tp = - (tp+ _ At/ (At + — At ) 

is tried. If this fails, the mid-point 
ip* = i p_+ (4> + - ip _)/2 

between 4 >_ and 4> + is tried. If this fails, the iterations are terminated 
and the coordinates are computed since a better value for >p is not 
obtainable. 


4. 5 THE COMPUTATION OF COORDINATES AND ACCELERATIONS 

The final values of the radius r and the function g have been 
obtained in the last iteration of Kepler's equation above. The functions 
S 1 and S 2 from the last iteration are used to determine 

(f-D = - «s 2 /r 0 
f = - <V<rr 0 ) 

(g -1) = - p s^/r . 
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Then the coordinates are obtained from the relations 


X = x 0 + [(f- 1) X Q + g ] 

y = y Q + t (f-i) y Q + g y Q ] 

z = Z Q + [ (f - 1) Z Q + g Z Q ] 
X = [f x Q + (g-1) X Q ] + x 0 

y = If y Q + (i- 1 ) y 0 J + 

z = [f Z Q + (g-1) z 0 l + z 0 


which give maximum accuracy with small values of r for which (f-1), 
g, f, (g-1) are each close to zero. Then 

3 


x 


Hx/r~ 


V = - ^y/r 

/ 3 

z = - fJL z/r 

give the acceleration components at time t, and the acceleration 
components 

*o = -f x o /r o 3 

v„ ■ 

•* I* j 3 

z o = z o /r o 

at time t are also computed, 

4, 6 THE COMPUTATION OF PARTIAL DERIVATIVES 


The partial derivatives are now computed by determining 

U = s 2 (t-t Q ) + M (c 4 - 3c 5 ) 4> 5 

and using locations in the matrix of partial derivatives to store 
elements of the matrices on the right side of the formulas 
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<— —v 

9 x/ 9 fj. 


r * n 

X X 


- S 2 /r 0 

9 y/ 

= 

y y 


U/r - s-. 

9 z/ 0^ 


z z 


L o 3-» 


— —> 

9 


1 • • * ~' 

XXX 


• -y< r r o r 

9 y/ 9ft 

= 

• If 

y y 


s /r 

9 z/9;u 


z z z 


2 0 


L- —' 


U/r - s„ 



L 0 3 J 


r— —i 

9 x~/dn 


• -* 

x o x o 


’ S 2 /r 

9 y Q /9n 

= 

y o y o 


-U/r + s 

9 z /9/i 


Z z 


L. 3 -J 

L 0 -J 


L o oj 



9 x q /9/h 


" x o 

x o 

x o 

9 y Q /3M 

= 

y o 

y o 

« # 

y o 

_ 9 Zq/9/^ 


_ z o 

z o 

*# 

z 

0 


- 

— 


V (r V 


S 2 /r 

- 

-u /r + s 3 ^ 


which give the partial derivatives with respect to ^ . Similarly, 
locations in the matrix of partial derivatives are used to store 
elements of the first matrix on the right side of the formulas for 
the 2x3 matrices 
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which are stored in locatiais of the inverse partial derivative matrix 
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Then the computations 


9x/ 3x 0 

3x/ 

9x/ 


— 

X 

x 


X 

1 -1 

:n° 

o 

:>s 

o 

tifj 

(f-l)+l 0 0 

a y/ a x Q 

9y/ 8y Q 

3y/& 0 


y 

» 

y 

h] + u 

y 

+ 

O 

I 

u-i 

O 

9 z / 9 X ^ 

d Z /dy Q 

9z / %J 


z 

z 


z 

L_ — 


0 0 (f-l) 

L +i - 


ax/9i Q 9x/9y o 9x/8z " 


X X 


X 

to *0 

~g 0 0 

3y/9x Q dy/ 9y Q 9y/9z Q 

= 

y y 

H - u 

y 

+ 

0 g 0 

9 z /9x 9z/9y„ 9z/9z 

0 0 0 


_ z z _ 


_z_ 


-1 

00 

o 

°l 


9x/ 9 x q 9x/ 9y Q 

9x/ 


X X 


X 

& 0 y 0 ^ 

f 0 0 

9y/9x Q 9y/9y Q 

9 z/9x^ 9 z/ 9y^ 

9y/9 z Q 

9z/ 3z 

oJ 


y y 

z z 

H + u 

** 

y 

z 

— 

+ 

0 f o 

_0 0 f 


9 x/ 9x^ 

9 x/9y Q 

9x/ 9z 

0 


— 

X 

X 

9y/ 9x q 

9y / 9y Q 

9y / 9z q 

= 

y 

0 

y 

9 z/ a^ 0 

9z/ 9y Q 

9 z/ 9 z q __ 


z 

» 

z 


give the 36 partial derivatives of x 



y> 


z, 



x, y, Z, with 



+ 


(g-l)+l 0 0 

o (g-l)+l 0 

0 0 (g-l)+l 


respect 


V V Z Q’ X Q’ V V Finall y* t ^ ie formulas 
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9 x Q / ck 

9x q / 9y 

ax Q / az 

> /- 

9 y Q /ax 

9y Q /9y 

9y Q / 9z 

- 

3 z;^/ 3 x 

9 z q/ ay 

9Zq/ 9 z 

- 

3 x^/ 3x 

9X Q / 9y 

9 x Q / 3 z 

r 

9 y Q /9x 

9 y Q /ay 

9y Q / 3z 

= - 

_3 z^/3x 

dz Q /dy 

9z q /9z _ 


3 Xq/3x 

ax Q / ay 

9x q / 3z 

' 

9 y Q /9x 

9 y 0 /ay 

9 y 0 / 9z 

- - 

^ Zq/3x 

di Q /dy 

dz^/d z 

- 

3 x^/3x 

dk Q /dy 

9x q /az 

< 

9 y Q / 9x 

dyj 9y 

dy Q /dk 

= 

3 z^/ 3x 

az Q /3y 

9z Q /3 z 



a x/3x o ax/ay Q 


9y/9x 0 9y/9y Q 


9z/9 x q dz/dy 


'9x/a: 


9z/ 3x 


0 


9x/ax 
9y / 9x 


0 


9y/9x 
8s/ 9x 


0 


ax/ay 


9y/ 9x 9y/ 9y 


dz/dy 


9x/3z r 


dy/di 


dz /az 


oj 


ax/ az 
9y/ 9z 


-> T 


0 


9z/az 


0- 


9i/ay. o 

9 x/9 z^ 

ay/ 9y 0 

ay / 9 zq 

di/dy 0 

9z/ 3z 

0 —1 

9x/ay Q 

9x/ az Q 

ay/ay Q 

9y/9 z q 

3z/9y Q 

3z/3 z 

oJ 

z o’ x o’ 

V z Q with 


0 

*ive the 36 inverse partial derivatives of x^, } 
respect to x, y, z, x, y, z. This concludes the computation of all 
quantities which are output by the subroutine. 
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4* 7 SOME PROPERTIES OF THE SUBROUTINE 


A Fortran 4 listing of the double-precision subroutine TWOBDY is 
given below. The compiled subroutine occupies less than 1500 single¬ 
precision storage locations on the IBM 7094. The minimum running 
time is a little over 10 milliseconds and occurs for the trivial case 
where the input value of r: is zero. The program will almost always 
run in less than 100 milliseconds for a non-zero value of ~T where 
no input approximation for ip is available. A good input approximation 
for $ can reduce the running time to as low as 15 milliseconds. Ex¬ 
tremely unusual cases with no input approximation for ip can take as 
much as 1000 milliseconds, but the solution will always be obtained. 

Accuracy checks indicate that errors in the outputs are always 
close to a lower bound below which accuracy improvements are un¬ 
warranted. The lower bound occurs because of uncertainties in the 
inputs due to the floating-point representation of numbers. For example, 
the absolute error in x should not be reduced below about 


| 3 x | 

• i 1016 * 0 1 + 

i _ 

ax 

1 

|l0- 16 y 0 M 

9x 

9x 0 

1 

9y o 

1 • 

9z o 

1 1 

. |io- 16 i 0 UI 

a 

X 

1 1 

‘OA.U 

l a x 

1 a io l 

a 

y o 

1 • 1 

1 9 z 

1 

i ax 

1. 1 io- 'Vi + 

1 • 

1 

|io' 

- 16 T [ 


1 

l x 

1 

T\ 


*ach of 

the eight inputs 

x , 
0 

V 

V V V 

V ' 


10-‘ 6 z 


, - 16 . 

10 z. 


known to only about sixteen significant figures. Partial derivatives have 
been checked by multiplyiig the 6x6 inverse partial derivative matrix. The 
difference between the product and the unit matrix is then easily ex¬ 
plainable in terms of errors caused by the floating-point representation 
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of the partial derivatives. 

Results obtained in some limited applications of the subroutine 
have been very satisfactory. The author will sincerely appreciate 
any suggestions, criticisms, or comments from those using the 
Fortran program. A Fortran version for the IBM 360 should be 
available shortly. There are several programs available which 
solve the equivalent problem accurately and efficiently for particular 
cases of two-body motion. The main advantage of the program de¬ 
scribed here is that it offers equal accuracy and efficiency with com¬ 
plete generality as well. This can become extremely important in 
applications where a compact program is desired to handle many 
different types of two-body motion. In any case, the generality is 
provided without any sacrifice of accuracy, efficiency, or compactness, 
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4. 8 THE FORTRAN 4 SYMBOLIC LISTING OF THE SUBROUTINE TWOBDY 


5IBFTC TWOBDY 

SUBROUTINE. TwOBDY(SOrTAUfMU»PSI*S»P»PIrPMU»P0MUrACC *ACCO»R#RO) 

GENERAL SOLUTION OF TWO BODY PROBLEM WITH PARTIAL DERIVATIVES 
FORTRAN 4 DOUBLE PRECISION SUBROUTINE FOR IBM 7094 WITH IBSYS SYSTEM 
SEE APRIL 1965 ASTRONOMICAL JOURNAL FOR FORMULATION BY W. H, GOODYEAR 

CALLING SEQUENCE IS AS FOLLOWS 

CALL TWOBDY <S0»TAUrMUrPSI»S»P»PI»PMU»P0MU»ACC»ACCO*R»R0) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TO 


DOUBLE PRECISION QUANTITIES IN CALLING SEQUENCE ARE AS FOLLOWS 
DOUBLE PRECISION SO(6)»TAU»MU*PSI 
l*S(6)»p(6»6)rPI(6f6> r PMU(6)»P0MU(6>*ACC(3)»ACCO(3)»R»R0 

INPUTS 

SO(l).SO(2),SO(3)=XOfYO*ZO=POSITION COMPONENTS AT REFERENCE TIME TO 
S0(4)»S0(5)»SO(6)=XDO*YDOrZDO=VELOCITY COMPONENTS AT REFERENCE TIME 

tau=time interval (T-to) from reference time to to solution time t 

MU=CONSTANT IN DIFFERENTIAL EQUATIONS (XDD»YnD»ZDD>=-MU*{X»Y»Z)/(R**3) 
PSI=APPR0XIMATI0N FOR FINAL SOLUTION PSI OF KEPLER’S EQUATION 

OUTPUTS 

PSI=GENERALIZED tCCENTRIC ANOMALY=SOLUTION OF KEPLERS EQUATION 
S(l)»S(2)»S(3)=X.Y»Z=P0SITI0N COMPONENTS AT SOLUTION TIME T=T0+TAU 
S(4)»S(5)*S(b)=Xu»YD»ZD=VELOCITY COMPONENTS AT SOLUTION TIME T=T0+TAU 
P(I.J)=PARTIAL DERIVATIVE DS(I)/DS0(J) OF S<I> WITH RESPECT TO S0(J) 

PI(I»J)=PARTIAL US0(I)/DS(J> WITH ROLES OF TO AND T REVERSED 

pmuiimpartial ds(i)/dmu of s(i) with respect to mu 

P0MU< IMPARTIAL DSO (1) /DMU WITH ROLES OF TO AND T REVEPSEO 
ACC<I)=-MU*S(I)/(R**3)=ACCELERATI0N COMPONENT AT SOLUTION TIME T 
ACCO(I)=-MU*S0(I)/(R0**3)^ACCELERATION COMPONENT AT REFERENCE TIME TO 
R=RADIUS AT TIME T=SQUARE ROOT OF<X**2+Y**2+Z**2) 

R0=RADIUS AT TIME T0=SQUARE ROOT OF(X0**2+Y0**2+Z0**2> 

ADDITIONAL DOUBLE PRECISION QUANTITIES FOR COMPUTATION 

2 »SIGO»ALPHA *PSIN»PSIP*A»AP»C0rCl»C2»C3»C4»C5X3>Sl»S2>S3»OTAU»DTAUN 
3»DTAUP#U» FMl» G r FD»GDMl 


C 

C START OF INITIAL COMPUTATIONS 

C COMPUTE RADIUS R0=SQUARE ROOT OF(X0**2+Y0**2+Z0**2) 
S1=0MAX1 (DABS (SOU) ) »DABS(S0(2) ) t DABS (SO (3)) ) 
S2=(SO(1)/Sl)**2+(SO(2)/Sl)**2+(S0(3)/Sl)**2 
R0=2. 


10 R=R0 
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R0=(R+S2/R)*.5 
IF(RO.LT.R) GO TO 10 
R0=R0*S1 

C COMPUTE OTHER PARAMETERS 

SIGO=SO(1)*SO<4)+SO(2>*SO(5)+SO(3)*SO(6) 
ALPHA=SO(4)**2+SO(5>**2+SO(6)**2-2.*MU/RO 
C INITIALIZE SERIES MOD COUNT M TO ZERO 
M=0 

C INITIALIZE BOUNDS PSIN AND PSIP FOR PSI OR SET PSI=0 IF TAU=0 
IF(TAU) 20»30»40 
20 PSIN=-1.0+38 
PSIP=0. 

DTAUN=PSIN 
DTAUP=-TAU 
GO TO 50 
30 PSI=0. 

GO TO 100 
40 PSIN=0. 

PSIP=+l.D+38 

DTAUN=-TAU 

DTAUP=PSIP 

C USE APPROXIMATION FOR PSI IF IT IS BETWEEN BOUNDS PSIN AND PSIP 
50 IF(PSI.GT.PSIN.AND.PSI.LT.PSIP) GO TO 100 
C TRY NEWTON’S METhOD FOR INITIAL PSI SET EQUAL To ZERO 
PSI=TAU/R0 

C SET PSI=TAU IF NEWTON’S METHOD FAILS 

IF(PSI.LE.PSIN.OR.PSI.GE.PSIP) PSI=TAU 
C END OF INITIAL COMPUTATIONS 
C 

C BEGINNING OF LOOP FOR SOLVING KEPLER’S EQUATION 
C BEGINNING OF SERIES SUMMATION 

C COMPUTE ARGUMENT A IN REDUCED SERIES OBTAINED BY FACTORING OUT PSI’S 
100 A=ALPHA*PSI*PSI 

IF(DABS(A).Le•1•) GO TO 120 

C SAVE A IN AP AND MOD A IF IT EXCEEDS UNITY IN MAGNITUDE 
AP=A 

110 M=M+1 
A=A*.25 

IF(DABS(A).61.1.) GO TO 110 
C SUM SERIES C5X3=3*S5/PSI**5 AND C4=S4/PSI**4 

120 C5X3=(l.+(l.+(l.+(l.-Ml.+(l.+(l.+A/342.)*A/272.)*A/2l0.)*A/15b.) 
1 *A/110.)*A/72.)*A/42.)/40. 

C4 =(1.+(l.+(1.♦(l.+(1.+(l.+(l.+A/306.)*A/240.)*A/182.)*A/132.) 
1 *A/90.)*A/56.)*A/30.>/24. 

C COMPUTE SERIES C3=S3/PSI**3»C2=S2/PSI**2»C1=S1/PSI»C0=S0 
C3=(.5+A*C5XB)/3. 
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C2= •b+A*C4 
Cl = 1•+ A*C3 
C0= 1•4A*C2 
IF(M.LE.O) GO TO 140 

C DEMOO SERIES CO AND Cl IF NECESSARY WITH DOUBLE ANGLE FORMULAS 
130 C1=C1*C0 

C0=2.*CU*C0-1. 

M=M-1 

IF(M.GT.O) GO TO 130 

C DETERMINE C2»C3'C4»C5X3 FROM C0»C1»AP IF DEMOD REQUIRED 
C2=(CO-1,)/Ap 
C3=(C1-1.)/AP 
C4=(C2-.b)/Af- 
C5X3=(3.*C3“.b)/AP 

C COMPUTE SERIES Sl»S2rS3 FROM C1»C2#C3 
140 S1=C1*PSI 

S2=C2*PS1*PS1 
S3=C3*PSI*PSI*PS1 
C END OF SERIES SUMMATION 

C COMPUTE RESIDUAL DTAU AND SLOPE R FOR KEPLER’S EQUATION 
G=R0*Sl+SIG0*S2 
DTAU=(G+MU*S3)-TAU 
R=DABS(R0*CO+(SIG0*S1+MU*S2)) 

IF(DTAU) 200 » 300»210 
C RESET BOUND 
200 PSIN=PSI 
DTAUN=DTAU 
GO TO 22u 
210 PSIP=PSI 
DTAUP=DTAU 

C TRY NEWTON'S METnOD AND INITIALIZE SELECTOR N 
220 PSI=PSI-DTAU/R 
N=0 

C ACCEPT PSI if IT IS BETWEEN BOUNDS PSIN AND PSIP 
230 IF(PSI.t>T.PS1N.AND.PSI.LT.PSIP) GO TO 100 
C SELECT ALTERNATE METHOD OF COMPUTING PSI OR STOP ITERATIONS 
N=N+1 

GO TO (1.2»3»4*300)»N 

C TRY INCREMENTING BOUND WITH DTAU NEAREST ZERO BY THE RATIO 4*DTAU/TAU 

1 IF(DABS(DTAUN).LT.DABS(DTAUP)) PSI=PSIN*(1.-(4.*DTAUN)/TAU) 
IF<DABS(dTAUP),LT.DABS(DTAUN)) PSI=PSIP*(1.-(4.*DTAUP)/TAU) 

GO TO 230 

C TRY DOUBLING BOUND CLOSEST TO ZERO 

2 IFiTAU.GT.O.) PSI=PSIN+PSIN 
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IFUAU.LT.O.) psi=psip+psip 

GO TO 250 

C TRY INTERPOLATION BETWEEN BOUNDS 

3 PSI=PSIN+(PSIP-PSIN)*(-DTAUN/(DTAUP-DTAUN)) 

GO TO 230 

C TRY HALVING BETWEEN BOUNDS 

4 PSI=PSIN+(PSIP-PSIN)*.5 
GO TO 230 

C END OF LOOP FOR SOLVING KEPLER’S EQUATION 
C 

C COMPUTE REMAINING THREE OF FOUR FUNCTIONS FMlfGfFDfGDM1 
300 FM1=-MU*S2/R0 
FD=-MU*S1/R0/R 
GDMl=-MU*S2/R 

C COMPUTE COORDINATES AT SOLUTION TIME T=T0+TAU 
DO 310 1=1 * 3 

S(I)=SQ(I)+(FM1*S0d) +G+S0 d+3)) 

S(I+3)=(FD*S0(I)+GDMl*S0(1+3))+S0(1+3) 

C COMPUTE ACCELERATIONS 

ACC(I)=-MU*S(I)/R/R/R 
310 ACCO(I)=-MU*S0(I)/RO/RQ/RO 
C END OF COMPUTATION FOR COORDINATES AND ACCELERATIONS 
C 

C COMPUTATION OF PARTIAL DERIVATIVES 
C COMPUTE COEFFICIENTS FOR STATE PARTIALS 

U= S2*TAU+MU*(C4-C5X3)*PSI*PSI*PSI*PSI*PSI 


Pdf 

1)=-(FD*S1+FM1/R0)/R0 

pdf 

2)=-FD*S2 

P(2f 

1)= FM1*S1/R0 

P (2 f 

2)= FM1*S2 

Pdf 

3)= P(lf2) 

pdf 

4)=-GDMl*S2 

P(2f 

3)= P(2f2) 

P(2» 

4)= G*S2 

P(3f 

1)=-FD*(C0/R0/R+1./R/R+l./RO/RO) 

P(3f 

2)=-(FD*Sl+GDMl/R)/R 

P(4f 

1)=-P(lfl) 

P(4f 

2)=-P(1f 2) 

P(3f 

3)= P(3f2) 

P(3f 

4)=-GDMl*Sl/R 

P (4 f 

3)=-P(lf2) 

P(4f 

4)=-P(1f4) 

C COMPUTE 

COEFFICIENTS FOR MU PARTIALS 

P(lf 

5)=-Sl/R0/R 

P(2f 

5)= S2/R0 

P(3f 

5)= U/R0-S3 
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P(l»6)=-P(l,b) 

P<2*6)= S2/R 
P(3»6)=-U/R+S3 
DO 4001=1,3 
C COMPUTE MU PARTIALS 

PMU(I)=-S<I> *P < 2.5)+S <I+3)*P(3.5) 

PMU(I + 3)= S( I)*P(1.5)+S(1+3)*P(2»5)+ACC (I)*P(3»5) 
P0MU<I)=-S0(i)*P(2.6)+S0(I+3)*P<3»6) 

P0MU(I+3)= S0(I)*P(1.6)+S0(I+3)*P(2»6)+ACC0(I)*P(3.6) 

C MATRIX ACCUMULATIONS FOR STATE PARTIALS 
DO 400 J-I.4 

PI(J»I>= P (J»1)*S0(I)+P(J*2)*S0(1+3) 

400 PI(J»I+3)= P(J»3)*S0(I)+P(J»4)*S0(1+3) 

DO 410 1=1»3 
DO 420 J=l»3 

P(I»J) =S(I)*PI(1»J) +S(I+3)*PI(2.J) +IJ*S< I+3)*ACC0(J) 

P(I»J+3) =S(I)*PI(1»J+3)+S(1 + 3)*PI(2»J+3)-U*S(1 + 3)*S0(J+3) 

P(I+3»j) =S( I) *PI (3. J) +S(I+3)*PI(4.J) +'J*ACC (I)+ACC0(J) 

420 P(1+3*J+3)=S(I)*PI(3.J+3)+S(1+3)*PI(4.J+3)-U*ACC(I)*S0(J+3) 
P(I»I) =P(I»I> +FM1+1. 

P(I * 1+3) —P(1*1+3) +G 

P<I+3*I) =P(I+3*I) +FD 

410 P(1+3*I+3)=P(I+3*I+3)+GDMl+l. 

C TRANSPOSITIONS FoR INVERSE STATE PARTIALS 
DO 430 1 = 1»3 
DO 430 J=1» 3 
PI(J+3.I+3)= P(I.J) 

PI(J+3.1) =-P(1 + 3 * J) 

PI(J.1+3) =-P(I.J+3) 

430 PI(J.I) = P(I+3.J+3) 

C END OF COMPUTATION FOR PARTIAL DERIVATIVES 
C 

C END OF PROGRAM - ALL OUTPUTS HAVE BEEN COMPUTEO 
RETURN 
END 


96 



REFERENCES 


Bower, E. , 1932, "Some Formulas and Tables Relating to Orbit 

Computation and Numeric Integration", Lick Observatory Bulletin 
No. 445, Vol. 16, p. 39. 

Goodyear, W. H. , 1965, "Completely General Closed-Form Solution 

for Coordinates and Partial Derivatives of the Two-Body Problem, " 
Astronomical Journal, Vol. 70, pp. 189-192. 

Herget,P., 1948, The Computation of Orbits, published privately by 
the Author, p. 48 and pp. 74-76. 

Herrick, S. I960, Astrodynamical Report No. 7, ASTIA Document 
AD 250 757, Par. 6Q. 

Kuhnert, F. , 1879, "Folgerungen aus v. Oppolzer's neuer Methode 

fur die Bearbeitung sp&terer Oppositionen, 11 Astronomische Nachrichten, 

Vol. 95, pp 145-149. 

Sconzo, P. , 1963, "Explicit Expressions for the 36 terms of the 

Jacobian Matrix Used in Orbit Computation," Mem. Soc. Astron. 

Halina, Vol. 34, No. 2. 

Sperling, H. , 1961, On the Computation of Keplerian Conic Sections, 

Report No. MTP-Aero-6l-19, Marshall Space Flight Center, Huntsville, 
Alabama. 

Stumpff,K. , 1947, "Neue Formeln und Hilfstafeln zur Ephemeridenrechnung, 

Astronomische Nachrichten, Vol. 275, p. 108. 

_, 1959, Himmelsmechanik (VEB Verlag, Berlin), Chapter V. 

, 1962, "Calculation of Ephemerides from Initial Values", 

NASA Technical Note D-1415. 

Sundmann, K. , 1912, "Memoire sur le probleme des trois corps", 

Acta Mathematika, Vol. 36, p. 127. 


NASA-Langley, 1966 CE-5^2 


97 



